SUMMARY

In this study, a straightforward and rapid methodology is proposed and tested to determine the seismic moment, the earthquake rupture length/duration and the static stress drop. To this purpose, three ground motion parameters, that is, P-wave peak acceleration (⁠|${P_a}$|⁠), velocity (⁠|${P_v}$|⁠) and displacement (⁠|${P_d}$|⁠) are evaluated as a function of time from the first P arrival. The average of the logarithm of the P-wave amplitude (LPDT curves), corrected for the distance-attenuation effect, is calculated using all the available stations in expanded P-wave time windows. The LPDT curves show an exponential growth shape and increase with time until they reach a constant value (plateau), which is related to the magnitude of the earthquake. From the obtained observations, we demonstrate that the corner time of the plateau level on the weighted-fit curve to the LPDT curves is related to the half-duration of the rupture. Thus, using the theoretical scaling, the source radius and stress drop can be obtained from the measured half-duration of the source. This method has been applied and tested to the records of the 2016–2017 Central Italy seismic sequence, with moment magnitude ranging between 3.4 and 6.5. Our study shows that source parameters match a self-similar, constant-stress-drop scaling with a relatively low average stress drop of about |$1.1 \pm 0.5\ \mathrm{ MPa}$|⁠, except for the largest event of the sequence showing a relatively higher stress release, which is associated with the dominant radiation from a localized high slip patch on the fracture surface. The proposed approach based on a simple time domain signal analysis is innovative and may complement longer spectral technique for fast estimating earthquake source properties.

INTRODUCTION

Undoubtedly, an earthquake is one of the most complex natural phenomena for which the accurate determination of all characteristic parameters is a difficult task. Although nowadays, estimating some earthquake properties such as location, magnitude and moment tensor is well-performed, characterizing the source properties such as average slip, fault length/surface and average stress drop is still a challenging issue (Allmann & Shearer 2009; Zollo et al. 2014; Kaneko & Shearer 2015).

The displacement spectra are generally modelled to obtain the average kinematic and static source parameters (Brune 1970; Madariaga 1976). To model the dynamics of the rupture and the related parameters, such as the dynamic stress drop, the high-frequency content of the spectrum is needed, and this may be affected by some complexities of the medium, such as distance attenuation and path effects (Allmann & Shearer 2009). Moreover, to simplify the modelling of the rupture process, the point-source approximation is usually assumed as a proxy for an extended fault surface. This is an acceptable assumption for small-to-moderate events (⁠|$M < 4.0\hbox{--}5.0$|⁠), while it may produce significant bias in the ground shaking prediction for large-magnitude events (⁠|$M \ge 6.0\hbox{--}6.5$|⁠) (Rydelek et al. 2007; Zollo et al. 2007; Festa et al. 2008).

Furthermore, the ground shaking prediction is generally based on the magnitude estimate and on the use of standard isotropic Ground Motion Prediction Equations that account only for the source-to-site distance and do not consider the azimuthal variation of the radiated wavefield. Thus, for large seismic events, the distance of the site of interest from a fault, the approximation of a single radiating point, may be inadequate. Nevertheless, there is strong evidence for a dominant finite-rupture effect on the distance/azimuth distribution of the ground shaking radiated by extended faulting phenomena especially at near-fault distances, for example, distance comparable with the fault length (Archuleta & Hartzell 1981; Somerville et al. 1997; Koketsu et al. 2016). Therefore, in the shake map computation, the finite extent of the source (length and width of the fault plane) is a relevant piece of information to be known, along with the focal mechanism that gives the fault geometry and orientation. Indeed, the knowledge of the expected rupture length and orientation will improve the accuracy in predicting the ground motion amplitude during moderate-to-large earthquakes, thus, producing realistic strong motion shake maps.

In the framework of EU SERA project (Seismology and Earthquake Engineering Research Infrastructure Alliance for Europe), different methodologies are being developed and tested to generate evolutionary ground shaking maps by considering a rupture kinematic description and reliable finite-fault model. In this regard, here we focus on the fault-size determination, through the refinement and testing of the methodology recently introduced by Colombelli & Zollo (2015), to rapidly estimate the event magnitude and the source parameters. In the original work, the authors computed the mentioned parameters (seismic moment and source duration) following the time evolution of the logarithm of P-wave peak displacement (⁠|${P_d}$|⁠) (called LPDT curve) using a data set of moderate-to-large Japanese events, with |${M_\mathrm{ w}}$| between 4.0 and 9.0.

In this study, following a similar approach, the source parameters (moment magnitude/average slip, fault length/surface and average stress drop) are automatically and rapidly determined for the events of 2016–2017 Central Italy seismic sequence (Chiaraluce et al. 2017). Here the original methodology is essentially modified to: (1) account for two other parameters (velocity and acceleration peaks, i.e. Pv and Pa) and (2) use a parametric approach where the observed curves are modelled with a three-parameter time function. The earthquakes occurred in an active seismic region (Meletti et al. 2016) of central Apennines in Italy and as it is shown in Fig. 1, almost all the events show an NW–SE striking, normal faulting mechanism. This long-lasting seismic sequence (Luzi et al. 2017) has started with the Amatrice earthquake, Mw 6.0 on August 24th while the largest event of the sequence, the Mw 6.5 Norcia earthquake, occurred on October 30th.

Data. The map shows the distribution of all the events of the 2016 Central Italy sequence (black open circles) and the epicentre position of the selected earthquakes (grey filled circles), with a variable size, depending on the magnitude. Dark grey triangles show the position of the stations used for the analysis. The focal mechanism solution (as provided by INGV) is also shown for the largest events (M ≥ 5.5).
Figure 1.

Data. The map shows the distribution of all the events of the 2016 Central Italy sequence (black open circles) and the epicentre position of the selected earthquakes (grey filled circles), with a variable size, depending on the magnitude. Dark grey triangles show the position of the stations used for the analysis. The focal mechanism solution (as provided by INGV) is also shown for the largest events (M ≥ 5.5).

While the approach proposed in this work has been originally conceived for the off-line characterization of the earthquake source properties, it is worth to note that the same approach could be also adapted to Earthquake Early Warning Systems and used for the real-time estimation of the main source parameters (Colombelli et al. 2012). The method is based on a time domain signal analysis and it is alternative or may complement more standard spectral techniques to estimate the seismic moment and source size/duration of earthquakes.

DATA AND METHODOLOGY

Database selection

We use a selection of earthquakes belonging to the 2016–2017 Central Italy seismic sequence, with moment magnitude ranging between 3.4 and 6.5. After a preliminary evaluation of all the available earthquakes, the original data set (135 earthquakes) is reduced to 28 events recorded by a minimum of 40 stations in the distance range between 0 and 100 km, including 12 events with moment magnitude larger than 4.7. We use a total number of 1895 of vertical components of the ground motion waveforms, recorded within the considered epicentral distance range. The map of the selected epicentres and stations is shown in Fig. 1. The selected stations belong to the Italian Strong Motion Network (Rete Accelerometrica Nazionale), operated by the Italian Department of Civil Protection, and to the Italian National Seismic Network, operated by the Istituto Nazionale di Geofisica e Vulcanologia (INGV). The P-phase arrival times have been manually picked on all the vertical acceleration waveforms and a 0.075 Hz high-pass Butterworth filter is applied to the displacement records to remove possible base-lines arising from low-frequency noise amplification due to two-times integration of the accelerometric records.

Parameter measurement and computation of LPDT curves

For each available waveform, we measure three initial peak amplitudes as the absolute maximum on acceleration, velocity and displacement waveforms, named |${P_a}$|⁠, |${P_v}$| and |${P_d}$|⁠, respectively. The peak amplitude parameters are measured in progressively expanded time window, starting from the P-wave onset using a time step of 0.01 s, which is twice the data sampling time (Fig. 2a). To correct the observed amplitude for the distance attenuation (Zollo et al. 2006; Nazeri et al. 2017), we first derived the empirical attenuation relationships for each parameter assuming the following polynomial form:
(1)
where M and R are the moment magnitude and hypocentral distance in km and |${P_x}$| is the considered ground motion quantity, with the subscript |$x$| denoting |$a$|⁠, |$v$| and |$d$| for acceleration in |$\mathrm{cm}\, {\mathrm{s}^{-2}}$|⁠, velocity in |$\mathrm{cm}\, {\mathrm{s}^{-1}}$| and displacement in |$\mathrm{ cm}$|⁠. To estimate the coefficients (⁠|${A_x}$|⁠, |${B_x}$|⁠, and |${C_x}$|⁠) of this relationship listed in Table 1, a least squared multiple regression analysis is used. We first estimate these coefficients for different time windows and find that the three coefficients do not change significantly with the window length, so that the same coefficients measured in a fixed time window of 3.0 s are used for the entire duration of the LPDT curves. Once the observed amplitudes are corrected (i.e. normalized to a reference distance of 1 km) by the hypocentral distance (hereinafter named |$P_a^C$|⁠, |$P_v^C$| and |$P_d^C$|⁠), using the equation |$\log \ P_x^C = \log {P_x} - \ {B_x}\log R$|⁠, the LPDT curves (in acceleration, velocity and displacement) are built as the average value of the logarithm of corrected peak amplitude parameters.
Relationships between LPDT curves and MRF. (a) The figure shows the different time windows (dashed lines) after the P-wave triggering on some sample traces to compute the peak amplitudes on each time window. (b) An example of final LPDT curve for a single event calculated from observed data. To model the observed curves, we used an exponential fitting function (expression 3), for which y0 (i.e. the intercept of the curve) is fixed to the first point of the curve. To correctly estimate ${P_L}\ $and ${T_L}$ we applied a weighted regression procedure to fit the LPDT curve (dotted curve). (c) The triangle-like moment rate function (MRF) and its relevant LPDT curve. The dashed line represents that the maximum level of the LPDT curves (Plateau time) occurs at the peak of the MRF (Peak Time).
Figure 2.

Relationships between LPDT curves and MRF. (a) The figure shows the different time windows (dashed lines) after the P-wave triggering on some sample traces to compute the peak amplitudes on each time window. (b) An example of final LPDT curve for a single event calculated from observed data. To model the observed curves, we used an exponential fitting function (expression 3), for which y0 (i.e. the intercept of the curve) is fixed to the first point of the curve. To correctly estimate |${P_L}\ $|and |${T_L}$| we applied a weighted regression procedure to fit the LPDT curve (dotted curve). (c) The triangle-like moment rate function (MRF) and its relevant LPDT curve. The dashed line represents that the maximum level of the LPDT curves (Plateau time) occurs at the peak of the MRF (Peak Time).

Table 1.

Coefficients of eq. (1) for each waveform.

ABCStandard error
Acceleration0.6−2.71.64± 0.37
Velocity0.76−2.3−1.26± 0.32
Displacement0.92−1.85−3.63± 0.35
ABCStandard error
Acceleration0.6−2.71.64± 0.37
Velocity0.76−2.3−1.26± 0.32
Displacement0.92−1.85−3.63± 0.35
Table 1.

Coefficients of eq. (1) for each waveform.

ABCStandard error
Acceleration0.6−2.71.64± 0.37
Velocity0.76−2.3−1.26± 0.32
Displacement0.92−1.85−3.63± 0.35
ABCStandard error
Acceleration0.6−2.71.64± 0.37
Velocity0.76−2.3−1.26± 0.32
Displacement0.92−1.85−3.63± 0.35
Table 2.

Summarizing the earthquake information and results from the presented study for the largest events (M ≥ 5.5).

DateTimeMagnitudeSource radius (km)Stress dropHalf-duration (s)
2016 August 241:36:326.006.500.763.14
2016 October 2619:18:055.907.171.703.47
2016 October 306:40:176.505.040.892.44
2017 January 1810:14:095.504.991.222.41
DateTimeMagnitudeSource radius (km)Stress dropHalf-duration (s)
2016 August 241:36:326.006.500.763.14
2016 October 2619:18:055.907.171.703.47
2016 October 306:40:176.505.040.892.44
2017 January 1810:14:095.504.991.222.41
Table 2.

Summarizing the earthquake information and results from the presented study for the largest events (M ≥ 5.5).

DateTimeMagnitudeSource radius (km)Stress dropHalf-duration (s)
2016 August 241:36:326.006.500.763.14
2016 October 2619:18:055.907.171.703.47
2016 October 306:40:176.505.040.892.44
2017 January 1810:14:095.504.991.222.41
DateTimeMagnitudeSource radius (km)Stress dropHalf-duration (s)
2016 August 241:36:326.006.500.763.14
2016 October 2619:18:055.907.171.703.47
2016 October 306:40:176.505.040.892.44
2017 January 1810:14:095.504.991.222.41

Following the approach of Colombelli et al. (2014) and Colombelli & Zollo (2015), the LPDT curves are computed at any time step after the P-wave arrival time and before the expected arrival of the S wave, which is set by a preliminary estimated empirical relationship between the SP traveltime and the hypocentral distance (see Fig. 2a). Thus, at each time step, the vertical component waveforms, possibly contaminated by the S-wave arrivals, are automatically excluded. Hence, by increasing the time window, the closest stations are eliminated one by one, and the computation is finally stopped when a minimum number of stations is available (five stations). Fig. 2(b) represents a typical generated LPDT curve using |$P_a^C$| parameter for the earthquake with magnitude 6.

Modelling of LPDT curves

To interpret the observed shape of LPDT curves, we simulate a triangle-like moment rate function (MRF), following the formulation of Sato & Hirasawa (1973) for a constant-velocity, circular rupture and generate its corresponding LPDT curve by using a similar procedure as for real data (see Fig. 2c). As it is clear in Fig. 2(c), the increase of the MRF corresponds to an increase of the LPDT curve and the beginning of the plateau on the LPDT curve (Plateau Time) occurs at the peak of the MRF. Thus, the plateau level carries information on the maximum amplitude of the MRF and the corresponding saturation time is related to the half-duration of the triangular function.

Theoretically, for a single station, the far-field source kinematic models indicate that the Pd versus time curve reproduces the shape of the apparent source time function (Madariaga 2015). The basic assumption of our work is that, by averaging Pd among many stations distributed over azimuth and distance, the resulting function approximates the MRF. In this representation, in order to reproduce the time scale of the MRF (time from P-wave onset at the source), we need to rearrange the recorded waveforms with a common time axis origin. This is the reason why the computation of Pd starts at each station at the arrival time of the P wave (Fig. 2a) and each point along the LPDT curve is the average amplitude (corrected by the distance), computed using the same time window for all the available stations.

If the source is represented by a triangle-like MRF, the maximum level of the LPDT curves occurs at the peak of the MRF. Theoretically, for the circular rupture model of Sato & Hirasawa, propagating at a uniform velocity from the nucleation of the rupture to the border, the relation between the source radius and half-duration is independent on the specific dynamic rupture model. Therefore, by considering all azimuthal coverage around the fault, the average half-duration (HD) is obtained as (Aki & Richards 2002):
(2)
where |$a$| is the source radius, |${v_r}$| is rupture velocity and |${v_p}$| is the P-wave velocity. Note that the P-wave velocity is fixed to 6 km s−1 and the rupture velocity varies between 0.5|${v_s}$| and 0.9|${v_s}$| where |${v_s}$| is the S-wave velocity. The final source radius is the average among all these computed values. Given an estimate of the half-duration of the source, the above equation allows estimating the radius of a circular earthquake rupture, without complex procedures and waveform analysis.
To model the observed LPDT curves, we adopted a specific function which is able to reproduce the initial increase of LPDT curves and their final plateau level. The function has the following form:
(3)
whose parameters are |${P_L}$|⁠, the plateau level of the curve; |${y_0}$|⁠, the intercept of the plot with the y-axis; and |${T_L}$|⁠, the time-constant that controls how rapidly an exponential function grows to the 63 per cent of its maximum value, that is,|${P_L}$|⁠.

RESULTS

The LPDT curves are obtained as the average among many stations distributed, in principle, over azimuth distance. A relevant aspect for the successful application of the proposed methodology is the good azimuthal coverage of recording stations and an adequate number of available data. Indeed, the use of very few data, or even single station curves, may result in irregular shapes of the LPDT curves, with the appearance of several steps from the initial time to the final plateau (Supporting Information Fig. S1). As it is clear from the plot (Supporting Information Fig. S1), the shape of the curve is far from the smooth average curve. Moreover, Single station-LPDT curves can be strongly affected by source/propagation effects (directivity, site effects), so that they cannot be used for a reliable representation of the source.

When there is a notable azimuthal gap in the seismic network or the number of available data is very limited (few stations), the data deficiency around the source may strongly affects the middle part of the LPDT curves, resulting in irregular shapes, with the appearance of intermediate, small steps before the final plateau value. Therefore, to avoid unrealistic estimation of the fitting parameters, we apply a two-step regression analysis, first a standard unweighted regression analysis and then a weighted fitting procedure for the initial part of the curves. The first analysis is applied to estimate preliminary values of |${P_L}$| and |${T_L}$| parameters. The fitting process is then repeated using a weighted fit, in which a larger weight (100 times larger) is assigned to the initial part of the curve (from the beginning to |${T_L}$|⁠) and a smaller weight is given to the remaining part of the curve. The quality of the weighted-fit curve is shown in Fig. 2(b).

Note that the parameter |${T_L}$| is only used in the second step of the fit process. Figs 3(a)–(c) show the LPDT curves in acceleration, velocity and displacement while the misfit values (at each time) are shown in bottom-right.

LPDT curves. (a–c) Average-logarithm of $P^C_a$, $P^C_v$ and $P^C_d$ in terms of different time windows exactly after the P-wave onset. In all panels, the mis-weighted fit is also shown (bottom-right) at the initial and last parts for all LPDT curves, computed based on the different used signals, that is, (a) acceleration, (b) velocity and (c) displacement.
Figure 3.

LPDT curves. (a–c) Average-logarithm of |$P^C_a$|⁠, |$P^C_v$| and |$P^C_d$| in terms of different time windows exactly after the P-wave onset. In all panels, the mis-weighted fit is also shown (bottom-right) at the initial and last parts for all LPDT curves, computed based on the different used signals, that is, (a) acceleration, (b) velocity and (c) displacement.

Fig. 4(a) shows the scaling of |${P_L}$| as a function of magnitude of the events. Thus, immediately after computing |${P_L}$| and estimating the magnitude based on the input data, the seismic moment |${M_0}$| can be easily calculated from the empirical relationship of the Hanks & Kanamori (1979) moment-magnitude scale. For each of the analysed events, the source radius is finally computed by inverting eq. (2), and assuming that the corner time of the plateau on the weighted fit curve is equal to the half-duration of the source.

Scaling relationships versus magnitude. (a) The plot shows the plateau values of the LPDT curves for acceleration, velocity and displacement as a function of magnitude. The best-fit line is shown by a black solid line. The best-fit linear regression equation is also shown on the plot. (b) Scaling of the logarithm of the source radius as a function of magnitude. The lines represent the theoretical scaling, with constant static stress-drop values (0.1, 1 and 10 MPa). The secondary y-axis presents the half-duration of the source time function. The inset plot shows the estimated stress drop for the individual earthquakes in the selected data set. The average stress-drop value ($\Delta \sigma \cong $$1.1 \pm 0.5\, \mathrm{ MPa}$) is shown as a dashed line. The grey squares represent the values computed by Madariaga's (1977) formula, using the average rupture length as obtained by two models (Cheloni et al. 2017; Xu et al. 2017). In both panels, the estimated parameters related to the LPDT curves of acceleration, velocity and displacement are shown with blue, green and red circles, respectively.
Figure 4.

Scaling relationships versus magnitude. (a) The plot shows the plateau values of the LPDT curves for acceleration, velocity and displacement as a function of magnitude. The best-fit line is shown by a black solid line. The best-fit linear regression equation is also shown on the plot. (b) Scaling of the logarithm of the source radius as a function of magnitude. The lines represent the theoretical scaling, with constant static stress-drop values (0.1, 1 and 10 MPa). The secondary y-axis presents the half-duration of the source time function. The inset plot shows the estimated stress drop for the individual earthquakes in the selected data set. The average stress-drop value (⁠|$\Delta \sigma \cong $||$1.1 \pm 0.5\, \mathrm{ MPa}$|⁠) is shown as a dashed line. The grey squares represent the values computed by Madariaga's (1977) formula, using the average rupture length as obtained by two models (Cheloni et al. 2017; Xu et al. 2017). In both panels, the estimated parameters related to the LPDT curves of acceleration, velocity and displacement are shown with blue, green and red circles, respectively.

DISCUSSION AND CONCLUSION

In this study, we have applied an optimized version of the method proposed by Colombelli & Zollo (2015) to a data set of earthquakes (⁠|$3.4 \le {M_\mathrm{ w}} \le 6.5$|⁠) occurred during the 2016–2017, Central Italy seismic sequence, including accelerometric data recorded at epicentral distances smaller than 100 |$\mathrm{ km}$|⁠. The proposed methodology is based on the use of LPDT curves with the main objective of calculating the seismic moment and source duration. As compared to the common procedures to compute the earthquake source parameters, including spectral domain parametric and non-parametric techniques, our method is straightforward, accurate and fast. Indeed, in the standard, spectrum-based approaches, the seismic moment is calculated from the low-frequency part of the displacement spectra, while the source radius is estimated from the corner frequency (Allmann & Shearer 2009). Stress drop (⁠|$\Delta \sigma $|⁠) is derived from moment (⁠|${M_\mathrm{ o}}$|⁠) and source radius (⁠|$a$|⁠) estimation using the Keilis-Borok (1959) equation, |$\Delta \sigma = \frac{7}{{16}}\frac{{{M_\mathrm{ o}}}}{{{a^3}}}$|⁠. In our methodology, the source parameters are computed from a simple measurement of the initial P-peak amplitude versus time, by considering the azimuthal averaging and the geometrical attenuation correction, using empirical scaling relationships (eq. 1).

According to the obtained results and the comparison with source parameter estimates using the other methods, which will be discussed later, no significant differences in uncertainties on source parameter estimations are observed. Indeed, although some uncertainties related to the manual phase-picking and data/fit-processing are predictable for our algorithm that may affect the results, the uncertainties of the spectrum-based methods such as the ones related to the bias between the corner frequency and attenuation parameter, affecting the spectral shape are not negligible (Kaneko & Shearer 2015).

As it is clear from Figs 3(a)–(c), for all the analysed events and for the three ground motion quantities (⁠|$P_a^C$|⁠, |$P_v^C$| and |$P_d^C$|⁠), the LPDT curves have a similar shape and scaling, with small initial values and a final plateau value, that is generally higher and reached in a longer time for the larger magnitude events. Therefore, the mentioned method can be applied to acceleration, velocity and displacement waveforms, with the main advantage of being easily exportable and adaptable to any kind of seismic network and of not requiring complicated data processing. For the three analysed parameters, we found a different scaling coefficient with magnitude, with a higher coefficient for displacement (around 1) and a smaller coefficient for velocity and acceleration (0.8–0.7). This progressively decrease from displacement to velocity/acceleration is consistent with what is predicted from the theory, in the assumption of a triangular shape of the Source Time Function (STF).

A unique filter (high-pass Butterworth filter) is only applied to the displacement signals in order to remove base-line effects on the time-series, while acceleration and velocity are used as they are recorded (only the mean value and the linear trend are removed). This is a further advantage of the methodology, since acceleration, velocity and displacement provide a complementary image of the entire spectral content of the source.

In the original work of Colombelli & Zollo (2015), the LPDT curves are modelled using a piecewise linear function, which is a too simple model to describe the continuous time evolution of the ground motion amplitude. Here a continuous and parametric exponential function to model the curves is adopted with only three parameters (⁠|${T_L},\ {P_L},\ {y_0}$|⁠): the characteristic time, |${T_L}$|⁠, and the plateau level, |${P_L}$|⁠, both controlling the evolution and shape of the curves. The scaling of |${P_L}$| is consistent with what has been found for the Japanese data set (Colombelli & Zollo 2015). The plateau level, |${P_L}$|⁠, is linearly increasing (in logarithmic scale) with the magnitude of the event, so that the larger the magnitude the higher the |${P_L}$| value (Fig. 4a). Therefore, as soon as the LPDT curve saturates, the magnitude is computed using eq. (1), that is, the empirical attenuation relationship between peak amplitudes, magnitude and the hypocentral distance.

Fig. 4(b) represents the estimated source length and half-duration of the source time function for different events in the data set which is computed from the corner time of the plateau level on the weighted-fit curve to the LPDT. In addition, we compute the theoretical scaling of the source radius as a function of magnitude using different fixed values of the stress drop from |$0.1$| to |$10\, \mathrm{ MPa}$| shown with different lines in Fig. 4(b). As it is evident from Fig. 4(b), for the considered magnitude range from small to moderate, the estimated values of the source radius are compatible with the theoretical expected trends and show a consistent self-similar, constant stress-drop scaling with magnitude. Specifically, for a magnitude 4.1, we find an average radius of about |$1.1 \pm 0.17\,\mathrm{ km}$|⁠, while |$6.5 \pm 0.76\, \mathrm{ km}$| is found for a magnitude 6.0 event. According to the computed source radius following the LPDT curves of the acceleration, velocity and displacement data, the average stress drops are |$1.2 \pm 0.6$|⁠, |$1.2 \pm 0.5$| and |$1.0 \pm 0.4$||$\mathrm{ MPa}$|⁠, respectively. It is worthwhile to mention that Bindi et al. (2004) using aftershocks of the 1997 Umbria–Marche seismic sequence (⁠|$1.4\ \le \ {M_\mathrm{ L}} \le \ 4.5$|⁠) in the central Italy have also found a self-similar scaling of static stress drop with the average about |$2.0 \pm 1.0\, \mathrm{ MPa}$| by analysing the S-wave spectra and using a non-parametric inversion approach. Moreover, the average stress drop estimated in this study is comparable with the average value of |$2.6\, \mathrm{ MPa}$| obtained for the 2009 L'Aquila sequence (Pacor et al. 2016).

For the four major earthquakes of the sequence with moment magnitude above 5.5, we also checked an independent estimation of the stress-drop value, considering the rupture size provided by two different models, based on the joint inversions of Synthetic Aperture Radar (SAR), Interferometric SAR and Global Position System data (Cheloni et al. 2017; Xu et al. 2017). To estimate the static stress drop of near-rectangular ruptures, we used the Madariaga (1977) equation to retrieve the static stress drop:
(4)
where |$W$| is the fault width, |$\mu $| is the rigidity at the source and |$C$| is a constant which depends on the specific fault geometry and slip direction. Madariaga (1977) evaluated |$C$| for circular, rectangular and elliptical fault geometries and found that it ranges between |${C_{\mathrm{ min}}} = \frac{{16}}{{7\pi }}\ = \ 0.73$| for circular ruptures and |${C_{\mathrm{ max}}} = \frac{\pi }{2}\ = \ 1.58$| for very long and thin ruptures. Note that the Keilis-Borok (1959) and Madariaga (1977) formulae provide the same static stress-drop value for circular fault ruptures. In our estimation of stress-drop values for the larger magnitude events of the Central Italy sequence, we used an intermediate value, for example, |$C\ = \ 1$|⁠, given the evidence of a near-rectangular faulting surface.

We find that the static stress drops (shown in Fig. 4b with the grey squares) for these events are in a good agreement with the values estimated in this study (listed in Table 2), with an average value of about |$1.5 \pm 0.5\, \mathrm{ MPa}$|⁠, except for the |${M_\mathrm{ w}} = \ 6.5$|⁠, Norcia event. Indeed, for the largest event of the sequence, the estimated value of the source radius (⁠|$5.0 \pm 0.9\, \mathrm{ km}$|⁠) is smaller than the expected average value (about 14 km) from the scaling relationship inferred from smaller magnitude events, with a consequent apparent higher value of stress drop (between about 15 and 44 MPa, with a mean value of |$27.6 \pm 15\, \mathrm{ MPa}$|⁠) and shorter half-duration of the source (with an average of 2.3 s versus a predicted value > 3.0 s). The high stress-drop value for this event is also consistent with the other independent estimates obtained from the energy-based procedures (Bindi et al. 2017; Picozzi et al. 2017). A possible reason for the underestimated source size and high stress drop of the largest event could be the effect of the high-frequency radiation from a dominant slip-patch (up to |$ >2\, \mathrm{ m}$| of average slip) (Cheloni et al. 2017) spreading over a smaller area than the final rupture surface. Although different authors reported a total rupture surface of about 80 to 200 |$\mathrm{ km}^2$| (Xu et al. 2017; Cheloni et al. 2017; Chiaraluce et al. 2017), the dominant slip-patch which has radiated during this earthquake covers a surface of about 65 |$\mathrm{ km}^2$| (Cheloni et al. 2017), which approximatively corresponds to the estimated source area (78 |$\mathrm{ km}^2$|⁠) in this study. In this case, we conclude that our methodology is sensitive to the seismic radiation from the dominant slip release fault patch and may be not able to retrieve the secondary and more complex effects of the total source time function.

The proposed method represents a simple and automatic approach to quickly estimate the earthquake magnitude and the expected length of the rupture, solely based on the continuous measurement of the initial P-wave peak amplitude. Although the required parameters are obtained here in an off-line analysis, the same methodology can be used to the future implementation of real-time, earthquake shaking prediction or even early warning. For instance, the earthquake rupture moment, length and stress release as an output of this method can be considered as initial reference for a source model to be used for computing the synthetic seismograms and rapid strong ground motion scenarios for earthquake impact evaluation. For real-time applications, the only piece of information needed is a reliable estimation of the earthquake location, in order to properly account for the path attenuation effect and normalizing the observed amplitudes. In terms of real-time applications, further analyses are needed to simulate the continuous data streaming, accounting for the P-wave propagation through the seismic network and to evaluate the real-time performance of the methodology. Assuming a standard velocity model, we can theoretically compute the time at which the measurements of required parameters could be available at the network. Fig. 4(b) shows that the source parameters would be available less than 1.0 s after the first P-wave arrival time for |$M < 4.0$|⁠, less than 2.5 s for |$4.0 < M \le 5.0$| and less than 4.3 s for two large events in the data set, that is, magnitude 6 and 6.5.

On the other hand, since the earthquake magnitude and the source properties are related to the released seismic energy, combining the S and P waves will obviously improve the final estimations but simultaneously arising some data processing complexity. While the use of the S wave is strongly related to the existence of a dense strong motion network around the region of interest, the automatic detection of the S phase is clearly not as simple as the identification of the P phase. Because, as mentioned above, the LPDT curve is an average of the P-wave amplitude corrected by the distance attenuation effect, if we want to add the S wave, we should also compute its attenuation relationship to correct the S-wave amplitude. Further computational complexities will be also due to the automatic selection of the coefficient for different phases. Thus, by only using the P wave, the process is more straightforward, rapid and noteworthy simple without any kind of complexities related to either considering the S waves or computing the spectra of the waves as a routine way to compute the rupture properties.

SUPPORTING INFORMATION

Figure S1. Examples of LPDT curves for a single event, computed in different ways. From top to bottom, the plot shows the LPDT curves, |$P^C_a$|⁠, |$P^C_v$| and |$P^C_d$|⁠, for a sample event with magnitude 6 (black curve). Green curves represent the LPDT curve for a single station while the red curves are generated by considering few stations, that is, with azimuth less than 30°.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

ACKNOWLEDGEMENTS

This work has been funded by EU SERA project (project ID 730900).

REFERENCES

Aki
K.
,
Richards
P.G.
,
2002
.
Quantitative Seismology
,
University Science Books
.

Allmann
B.P.
,
Shearer
P.M.
,
2009
.
Global variations of stress drop for moderate to large earthquakes
,
J. geophys. Res.
,
114
,
B01310
,
doi:10.1029/2008JB005821
.

Archuleta
R.J.
,
Hartzell
S.H.
,
1981
.
Effects of fault finiteness on near source ground motion
,
Bull. seism. Soc. Am.
,
11
,
937
957
.

Bindi
D.
,
Castro
R.R.
,
Franceschina
G.
,
Luzi
L.
,
Pacor
F.
,
2004
.
The 1997–1998 Umbria-Marche sequence (Central Italy): source, path and site effects estimated from strong motion data recorded in the epicentral area
,
J. geophys. Res.
,
109
,
B04312
,
doi:10.1029/2003JB002857
.

Bindi
D.
,
Spallarossa
D.
,
Picozzi
M.
,
Scafidi
D.
,
Cotton
F.
,
2017
.
Impact of magnitude selection on aleatory variability associated with ground-motion prediction equations: Part I—Local, energy, and moment magnitude calibration and stress-drop variability in Central Italy
,
Bull. seism. Soc. Am.
,
108
,
1427
1442
.

Brune
J.
,
1970
.
Tectonic stress and spectra of seismic shear waves from earthquakes
,
J. geophys. Res.
,
75
,
4997
5009
.

Cheloni
D.
et al. .,
2017
.
Geodetic model of the 2016 Central Italy earthquake sequence inferred from InSAR and GPS data
,
Geophys. Res. Lett.
,
44
,
6778
6787
.

Chiaraluce
L.
et al. .,
2017
.
The 2016 Central Italy seismic sequence: a first look at the mainshocks, aftershocks, and source models
,
Seismol. Res. Lett.
,
88
,
757
771
.

Colombelli
S.
,
Zollo
A.
,
2015
.
Fast determination of earthquake magnitude and fault extent from real-time P-wave recordings
,
Geophys. J. Int.
,
502
,
1158
1163
.

Colombelli
S.
,
Zollo
A.
,
Festa
G.
,
Picozzi
M.
,
2014
.
Evidence for a difference in rupture initiation between small and large earthquakes
,
Nat. Commun.
,
5
,
3958
.

Festa
G.
,
Zollo
A.
,
Lancieri
M.
,
2008
.
Earthquake magnitude estimation from early radiated energy
,
Geophys. Res. Lett.
,
35
,
L22307
.

Hanks
T.C.
,
Kanamori
H.
,
1979
.
A moment magnitude scale
,
J. geophys. Res.
,
84
(
5
),
2348
2350
.

Kaneko
Y.
,
Shearer
P.M.
,
2015
.
Variability of seismic source spectra, estimated stress drop and radiated energy, derived from cohesive-zone models of symmetrical and asymmetrical circular and elliptical ruptures
,
J. geophys. Res.
,
120
,
1053
1079
.

Keilis-Borok
V.I.
,
1959
.
On estimation of the displacement in an earthquake source and of source dimensions
,
Ann. Geofis.
,12
,
205
214
., .

Koketsu
K.
et al. .,
2016
.
Widespread ground motion distribution caused by rupture directivity during the 2015 Gorkha, Nepal earthquake
,
Sci. Rep.
,
6
,
28536
.

Luzi
L.
et al. .,
2017
.
The Central Italy seismic sequence between August and December 2016: analysis of strong-motion observations
,
Seismol. Res. Lett.
,
88
(
5
),
1219
1231
.

Madariaga
R.
,
1976
.
Dynamics of an expanding circular fault
,
Bull. seism. Soc. Am.
,
66
(
3
),
639
666
.

Madariaga
R.
,
1977
.
Implications of stress-drop models of earthquakes for the inversion of stress drop from seismic observations
,
Pure appl. Geophys.
,
115
,
301
316
.

Madariaga
R.
,
2015
.
Seismic source theory
, in
Treatise on Geophysics
, 2nd edn, pp.
51
71
., ed.
Gerald
S.
,
Elsevier B.V

Meletti
C.
,
Visini
F.
,
D'Amico
V.
,
Rovida
A.
,
2016
.
Seismic hazard in central Italy and the 2016 Amatrice earthquake
,
Ann. Geophys.
,
59
,
doi:10.4401/AG-7248
.

Nazeri
S.
,
Shomali
Z.H.
,
Colombelli
S.
,
Elia
L.
,
Zollo
A.
,
2017
.
Magnitude estimation based on integrated amplitude and frequency content of the initial P wave in earthquake early warning applied to Tehran, Iran
,
Bull. seism. Soc. Am.
,
107
(
3
),
1432
1438
.

Pacor
F.
et al. .,
2016
.
Spectral models for ground motion prediction in the L'Aquila region (central Italy): evidence for stress-drop dependence on magnitude and depth
,
Geophys. J. Int.
,
204
(
2
),
697
718
.

Picozzi
M.
,
Bindi
D.
,
Brondi
P.
,
Di Giacomo
D.
,
Parolai
S.
,
Zollo
A.
,
2017
.
Rapid determination of P wave-based energy magnitude: in-sights on source parameter scaling of the 2016 central Italy earthquake sequence
,
Geophys. Res. Lett.
,
44
,
4036
4045
.

Rydelek
P.
,
Wu
C.
,
Horiuchi
S.
,
2007
.
Comment on “Earthquake magnitude estimation from peak amplitudes of very early seismic signals on strong motion records” by Aldo Zollo, Maria Lancieri, and Stefan Nielsen
,
Geophys. Res. Lett.
,
34
,
L20302
,
doi:10.1029/2007GL029387
.

Sato
T.
,
Hirasawa
T.
,
1973
.
Body wave spectra from propagating shear cracks
,
J. Phys. Earth
,
21
,
415
431
.

Somerville
P.G.
,
Smith
N.F.
,
Graves
R.W.
,
Abrahamson
N.A.
,
1997
.
Modification of empirical strong ground motion attenuation relations to include the amplitude and duration effects of rupture directivity
,
Seismol. Res. Lett.
,
68
,
199
222
.

Xu
G.
,
Xu
C.
,
Wen
Y.
,
Jiang
G.
,
2017
.
Source parameters of the 2016–2017 Central Italy earthquake sequence from the Sentinel-1, ALOS-2 and GPS data
,
Remote Sens.
,
9
,
1182
,
doi:10.3390/rs9111182
.

Zollo
A.
,
Lancieri
M.
,
Nielsen
S.
,
2006
.
Earthquake magnitude estimation from peak amplitudes of very early seismic signals on strong motion records
,
Geophys. Res. Lett.
,
33
,
L23312
,
doi:10.1029/2006GL027795
.

Zollo
A.
,
Lancieri
M.
,
Nielsen
S.
,
2007
.
Reply to comment by P. Rydelek et al. on ‘‘Earthquake magnitude estimation from peak amplitudes of very early seismic signals on strong motion records’’
,
Geophys. Res. Lett.
,
34
,
L20303
,
doi:10.1029/2007GL030560
.

Zollo
A.
,
Orefice
A.
,
Convertito
V.
,
2014
.
Source parameter scaling and radiation efficiency of microearthquakes along the Irpinia fault zone in southern Apennines, Italy
,
J. geophys. Res.
,
119
,
3256
3275
.

Colombelli
S.
,
Zollo
A.
,
Festa
G.
,
Kanamori
H.
,
2012
.
Early magnitude and potential damage zone estimates for the great Mw 9 Tohoku-Oki earthquake
,
Geophys. Res. Lett.
,
39
,
L22306
,
doi:10.1029/2012GL053923
.

This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://dbpia.nl.go.kr/journals/pages/open_access/funder_policies/chorus/standard_publication_model)

Supplementary data