A rapid response magnitude scale for timely assessment of the high frequency seismic radiation

In this work the scaling of seismic moment (M0) and radiated energy (Er) is investigated for almost 800 earthquakes of the 2016–17 Amatrice-Norcia sequences in Italy, ranging in moment magnitude (Mw) from 2.5 to 6.5. The analysis of the M0-to-Er scaling highlights a breaking of the source self-similarity, with higher stress drops for larger events. Our results show the limitation of using M0, and in turn Mw, to capture the variability of the high frequency ground motion. Since the observed seismicity does not agree with the assumptions on stress drop in the definition of Mw, we exploit the availability of both Er and M0 to modify the definition of Mw and introduce a rapid response magnitude (Mr), which accounts for the dynamic properties of rupture. The new Mr scale allows us to improve the prediction of the earthquake shaking potential, as shown by the reduction of the between-event residuals computed for the peak ground velocity. The procedure we propose is therefore a significant step towards a quick assessment of earthquakes damage potential and timely implementation of emergency plans.

In the aftermath of an earthquake, the rapid assessment of both location and extension of potentially damaging ground shaking is a primary task for seismological agencies supporting emergency managers. In this context, shaking maps 1 become a de-facto standard for a timely dissemination of the ground shaking experienced in the area struck by an earthquake. To improve the spatial resolution of such maps for rapid response actions, the rapid determination of an earthquake size and location supports the information provided by the actual ground motion measurements (if available) and predicted ones. The moment magnitude M w 2,3 is used by the seismological community as the primary measure of the earthquake size. Since M w is based on an estimate of the seismic moment M 0 4,5 , it provides fault-averaged, low-frequency information on source processes but relatively less information about the small-wavelength high-frequency rupture details 6 . For example, earthquakes with similar M w but different stress drop (Δσ) can generate different ground motion levels 7 , suggesting that a rapid assessment of Δσ could allow more reliable predictions of the earthquake-induced ground motion severity to engineering structures (hereinafter, shaking potential). Moreover, recent studies [8][9][10] showed that the Δσ variability is a key parameter for explaining the between-event residuals at short periods. Since different magnitude scales provide different information about the static and dynamic features of the earthquake rupture, magnitudes other than M w could better characterize the earthquake size in terms of high-frequency energy release [11][12][13] . For example, M w can be complemented with a magnitude scale based on the high-frequency level of the Fourier spectrum 14 . Using teleseismic broadband P-wave recordings, a magnitude scale (M e ) was introduced 15 based on measurements of the radiated energy (E r ) and revising the Gutenberg and Richter relationship 16 between E r and the surface-wave magnitude M s . Since M e is directly linked to the source dynamics, it is more sensitive to high-frequency source details such as variations of the slip and/or stress conditions, and the dynamic friction at the fault surface during the rupture process. Indeed, high energy-to-moment ratios indicate that the intensity of radiated energy at high frequencies is large relative to the size (measured by moment) of an earthquake, with significant implications on hazard assessment 17 . Although automatic procedures for the rapid estimation of M e using P-wave recordings have been proposed both at teleseismic 18 and local 19 distances, a strategy to combine the information provided by M w and M e for a rapid assessment of the damage potential of an earthquake has not been proposed yet.
In this study, we present a new procedure to measure the earthquake size using rapid assessments of both E r and M 0 , considering S-wave recordings within 100 km from the epicenter. The methodology is applied to Central Italy, used as an example region where the seismic hazard for residential buildings is dominated by close-distance earthquakes of low-to-moderate magnitude 20 (i.e., from 4.5 to 6.5). We first calibrate empirical attenuation relationships between the integral of squared velocity (IV2 S ) and the peak displacement (PD S ) with respect to E r and M 0 , respectively, considering 229 earthquakes ranging from M w 2.4 to 6.1, mostly belonging to the L' Aquila (2009) seismic sequence 21 . Then, the procedure is applied to about 775 earthquakes of the 2016-2017 Central Italy sequence 22 with M w in the range 2.5-6.5. Innovatively, we use the estimates of E r and M 0 (i.e., through the slowness parameter 23 ) to introduce a rapid response magnitude (M r ). The new M r magnitude scale is tied to the non-saturating M w , but it includes a term which accounts for the difference between the earthquake dynamic conditions and those assumed by Kanamori for the original definition of M w . Finally, the impact of the new magnitude scale is discussed in terms of assessment of the earthquake shaking potential quantified through the peak ground velocity (PGV) and peak ground acceleration (PGA).

Background
The moment magnitude M w is defined as 2 . . M (logM 9 1)/1 5 (logM where the slope value of 1.5 and the constant −4.8 are inherited from the relationship between log(E r ) and the surface-wave magnitude M s 16 , whilst the term −4.3 derives from the following assumptions: (a) the energy required for fracturing is negligible; (b) the final average stress and the average stress during faulting are equal (also known as "complete stress drop" Δσ or "Orowan's model" 24 ); (c) the average rigidity in the source area (μ) is ranging from 3 to 6 × 10 4 MPa under average crust-upper mantle conditions; (d) Δσ is nearly constant for very large earthquakes, with values between 2 and 6 MPa.
The above assumptions can be also formalized as: Introducing the slowness parameter term θ = log(E r /M 0 ) 23 , equation (2) corresponds to: K that is, θ assumes the value −4.3 under the Kanamori assumptions 2 on Δσ and μ. However, it is well-known that rupture processes may deviate from the complete stress drop, ranging between the "partial stress drop" 25 and "frictional overshoot" 26 models. Global datasets of earthquakes with M w varying between 5.5 and 9 suggest that θ varies between −4.7 and −4.9 15,27,28 , which would correspond to an average global stress drop a factor 3 to 4 smaller than that assumed by Kanamori 2 . This difference in the average stress drop values is also reflected by the expected over estimation of about 0.27 m.u. of M e with respect to M w when an earthquake satisfies the condition θ = θ K

Dataset.
In this study, we analyze about 1000 earthquakes with M w between 2.5 and 6.5, occurred in Central Italy between July 2008 and September 2017 (see Text Sa). This area of the Apennines is characterized by southwest-northeast extension (i.e., active NE-and SW-dipping normal and normal-oblique faults) since the Middle Pliocene 30 . The dataset includes the main seismic sequences that have occurred over the past 10 years in Central Italy: the 2009 L' Aquila sequence 21 , and the 2016-2017 Central Italy sequence 31 (Fig. 1a). The four largest analyzed The data were recorded by the Italian strong motion network (RAN) and by the permanent and temporary stations of other networks (see Text Sa and Fig. 1). To reproduce a rapid response procedure, the recordings are pre-processed and the squared-velocity parameter (IV2 S ) and the peak displacement (PD S ) are both measured on direct S-waves following an analysis scheme suitable for real-time operations 32 , (see Text Sb). From the whole compiled dataset, which consists of more than 1400 earthquakes, we extracted a subset considering the following selections: hypocentral distance smaller than 100 km; events recorded by a minimum of 8 stations and stations having at least 8 records; the sum of SNR for the three components ≥200. The selected dataset is composed by 1004 earthquakes recorded at 340 stations. This dataset is split in two parts. The first part consists of 229 earthquakes, which occurred between July 2008 and January 2016 (hereinafter referred to as dataset D1), including the 2009 L' Aquila M w 6.3 sequence. Data set D1 is used to calibrate the attenuation models used in this study. The second dataset (hereinafter referred to as D2) is composed of 775 earthquakes occurred in the study area since January 2016 and it is used to exemplify the application of the proposed methodologies. For D2, the number of good quality recordings per event varies from 10 to 50 ( Figure S1).

Results
Rapid response assessment of radiated energy and seismic moment. The rapid assessment of E r and M 0 is based on extracting, from the time histories, proxies for these two source parameters, and correcting them for attenuation along the path. In particular, the radiated seismic energy is estimated from the squared velocity integrated over the S-wave time window 33 (IV2 S ), whereas the seismic moment is derived from the S-wave peak-displacement 34 (PD S ). The parameters IV2 S and PD S are linked to E r and M 0 through empirical attenuation models derived using data set D1, as detailed in the section Method (eqs 8 and 9). Since linear non-parametric regressions 19 are applied, the attenuation models are tables listing the attenuation coefficients for the discretized distance range. The attenuation models are reported in Table S1.
The attenuation models (8) and (9) calibrated over data set D1 are used to compute E r and M 0 of earthquakes included in D2. The E r and M 0 values are obtained by averaging over several recording stations (see Text Sc and Figure S2 for details). Figure (2) shows that a linear scaling between E r and M 0 holds over 7 orders of magnitude for seismic moment. The observed M 0 -to-E r scaling well compares to that of other tectonic areas ( Figure S3). The best-fit model is:  Energy to moment scaling and event specific adjustment to M w . A straightforward application for rapid response purposes of the E r and M 0 estimates derived for the 2016-2017 dataset would consist of converting them into energy and moment magnitudes. In the case of M 0 , we adopt the original definition proposed by Kanamori 2 (i.e., Eq. 1). By contrast, in the case of the energy magnitude, we use the calibration dataset to parameterize a linear relationship between E r and the local magnitude (M L ) 35 . The best-fit model obtained by an iteratively reweighted least squares 36 and with 200 bootstrap replications 37 considered to assess the uncertainties ( Figure S4) is: L r with R 2 equal to 0.93 and the standard deviation of the residuals is equal to 0.13. Applying Eq. (5), the derived magnitudes agree with M L , but being tied to E r , can be extended toward larger earthquakes to avoid magnitude saturation effects 38 . The magnitudes obtained by the 2016-2017 E r estimates are referred to as energy-based local magnitude 19 (M le ). Figure (3a) shows that, while M w and M le well agree for M w > 5, for smaller events M w tends to progressively overestimate M le , and the (M w − M le ) difference increases with decreasing stress-drop (Δσ) 39 .
The M 0 -to-E r scaling (Fig. 2) and the differences between the energy and moment magnitudes in relation to Δσ (Fig. 3a) reflect a regional behavior of the source scaling which deviates from the Kanamori's assumptions θ K = −4.3 (Eq. 3). As shown in Fig. (2), the ratio log(E r /M 0 ) is expected to be smaller than the values implied by θ K = −4.3 for M w < 4.3. Considering the magnitude distribution of the analyzed earthquakes, most of the empirical θ values are smaller than θ K Fig. 3b, in agreement with Δσ smaller than the value assumed by Kanamori 2 . Therefore, these earthquakes can potentially generate a ground shaking variability at high frequency larger than the variability expected by assuming θ constant as in Eq. (1).
In the context of early warning and rapid response applications where event-specific analysis is of concern, it is worth considering a new magnitude scale that accounts for the differences in the events dynamic conditions with respect to those assumed by Kanamori for the original definition of M w . Figure (3c) shows that ΔM = (M le − M w ) scales with Δθ = (θ − θ K ) and correlate with Δσ While Δθ provides a mean for inferring, in the rapid response timeframe, how much the occurring earthquake deviates from the Kanamori's condition (i.e., if the earthquake is as energetic as it is expected from the θ K condition), an empirical linear model between Δθ and ΔM can be used for modifying the original definition of M w : Δ = . ± . Δθ − . ± . M (0 34 0 04) (0 48 0 04); ( 6) with R 2 equal to 0.88 and a standard deviation, assessed by a bootstrap approach 37 , equal to 0.16. Therefore, given Eq. (6), ΔM can be added to M w for deriving a new rapid response magnitude scale: where M r is used in the remainder of this work to indicate the new magnitude scale introduced to take into account the event-specific deviation of θ from θ K . Figure (3d) compares M r with M w estimated by applying a non-parametric spectral inversion approach 39 . It is worth noting that for earthquakes characterized by Δσ values from ~5 to ~15 MPa, ΔM varies between −0.2 to 0.2 magnitude units. Therefore, the difference between M r and M w is negligible for the largest events of the sequence (i.e. for M w > 5) which fulfill Kanamori's condition. However, for smaller earthquakes characterized by Δσ below 2 MPa, the application of the ΔM factor leads to M r being significantly smaller than M w (Fig. 3d).
To quantify the impact of M r on predicting the ground shaking potential, we follow the earthquake engineering approach of analyzing the event-specific deviations from the expected mean regional attenuation trend. We refer to PGV as a measure of earthquake shaking potential, since it is controlled by the seismic energy at intermediate frequencies and, hence, better suited than peak ground acceleration (PGA) for statistical analysis correlating ground motion with damage 40 . We calibrate a ground motion prediction equation (GMPE) to model the PGV scaling with distance and magnitude in the target area (see Text Sc). The calibration is repeated twice, considering either M w or M r . Then, the between-event residuals 41 (δBe) are computed as the average difference, for any given event, between the PGV measured at different stations and the corresponding values predicted by the GMPE. However, considering that PGA can be anyway of interest for engineering 42 or seismological 10,43,44 applications, the same analyses are repeated for PGA, (see Text Sc). Figure (4a) shows that the standard deviation of δBe computed considering PGV reduces from 0.18 to 0.08 when M w is replaced with M r , meaning that M r better accounts for those variations in the rupture process which can introduce systematic event-dependent deviations from the mean regional PGV scaling. As shown in Fig.  (4b), the reduced δBe variability is mostly related to the small magnitude events (i.e. M w < 4.5) which show the largest deviation from the θ = θ K condition. A similar result is obtained also for PGA ( Fig. 4c and d), for which the standard deviation of δBe reduces from 0.22 to 0.13 when M w is replaced by M r .

Discussion
This study was motivated by the need of developing a procedure for the rapid assessment of the earthquake shaking potential. To accomplish such a task, we worked on two topics: first, we introduced an approach for the rapid assessment of the seismic moment and the radiated energy; second, we used these two source parameters to define a measure of the earthquake size which is not only informative for the average tectonic effects but it accounts for the efficiency of the earthquake to radiate seismic energy at high-frequency.
Concerning the first topic, standard approaches adopted by seismic observatories for estimating M 0 from moment tensor analysis 45,46 are delayed by the requirement of acquiring full waveforms at regional-teleseismic distances. For instance, in Italy the moment tensor solutions are provided by INGV generally between 6 to 9 minutes after the earthquake occurrence (http://info.terremoti.ingv.it/en/help#TDMT) and ShakeMaps are released to the Civil Protection after this time span (http://shakemap.rm.ingv.it/shake/about.html). In this work, we have proposed an approach for quick and robust estimations of M 0 and E r based on proxies measured on S-waves recorded by a local network. The application of the attenuation models calibrated for both M 0 and E r to 775 earthquakes of the 2016-2017 Central Italy sequence confirmed the suitability of the approach for rapid-response applications. Figure (S8) shows the comparison between the E r estimates obtained by the procedure proposed in this study, which follows Picozzi et al. 19 , and the E r obtained for the same dataset by a non-parametric inversion approach 39 , which follows Kanamori et al. 38 . As shown within Figure (S8), the two considered procedures provide consistent results, within a factor 1.5/2 (i.e., the log(E r ) differences are within 0.2/0.3 log 10 units). We observe that the rapid estimates of E r tend to be slightly larger than those computed by Bindi et al. 39 . At least part of this overestimation could be ascribed to local site amplification effects, which have been considered in Bindi et al. 39 through site correction terms, not included in equation (8), an issue of interest for future studies. However, considering that E r estimates on the same earthquakes provided by different investigators often differ up to one order of magnitude 47 , the result supports the reliability of the rapid estimates of E r from IV2 S . The availability of a dense network in the area under study was certainly an important condition to average out radiation pattern and rupture directivity effects on E r and M 0 estimates. Future examination will explore the bias on E r and M 0 estimates for network configurations less dense than the one considered.
Regarding the earthquake damage potential assessment, the analyzed data set represented a suitable test-case to discuss the limitation of using M 0 to capture the variability of the high frequency ground motion. The topic is so important that it has been recently debated by different studies [7][8][9][10] . In this study, we have considered the same scientific topic but from the point of view of an observatory which must provide real-time information to support civil protection decisions and seismic risk mitigation actions. It is important to note that early-warning and rapid-response systems often deal with seismic sequences occurring in a specific target region. Since the within-sequence stress drop variability is expected to be larger than the between-mainshock stress drop variability 48  the earthquakes analyzed in this study varies from about 0.1 to 15 MPa, leading to slowness ratios significantly different from θ K (Fig. 3b).
The rapid assessment of M 0 and E r can thus be exploited to define a new rapid response magnitude scale (M r ) tied to M w but including a term that accounts for the difference between the actual energy-to-moment ratio and the value used by Kanamori 2 .
To demonstrate the potential value of the rapid response magnitude M r in better capturing the event-to-event shaking potential variability, we quantified the reduction of the between-event standard deviation for PGV and PGA, which is controlled by the Δσ variability. Considering Δσ values 37 and the scheme proposed by Choy 49 over the energy-to-moment ratio, we observe that earthquakes with M w > 4 (showing Δσ from 5 to 15 MPa) are moderately or highly energetic, whilst small magnitude events (characterized by Δσ < 1 MPa) are enervated ( Figure S5). In our opinion, such large variability in the source properties justifies the introduction of the new magnitude scale for rapid response procedures, which is capable of better predicting the ground motion over a frequency range of engineering interest.
The proposed approach for computing E r , M 0 and M r can be implemented into algorithms for real-time operations 50,51 , so that M r estimates can be available as soon as the S-waves are recorded by a minimum number of stations that allow stable E r and M 0 estimates to be obtained (i.e., within a few tens of seconds, depending on the seismic network density and telemetry efficiency), and can be exported to other areas for which the attenuation model for E r and M 0 can be calibrated. The application of the procedure for magnitudes larger than those in our dataset might be influenced by the saturation limit of the proxies used for the estimation of E r and M 0 . We expect that extending the observation scale from local to regional networks and considering other seismic phases than direct S-waves would allow for estimates of M r to be extended to earthquakes with magnitude larger than M w 6.5. The saturation limit for the rapid-response procedure here proposed will be investigated in future studies.
The implications of our results are for a broad scientific audience (e.g., seismological, engineering seismology, disaster manager communities). Indeed, the implementation of our procedure for computing the seismic moment (M 0 ), the radiated energy (E r ) and M r in real-time operations can allow civil protection operators to  better assess within a rapid response timeframe (i.e., within a few tens of seconds from the earthquake's origin time) an earthquake damage potential and the timely implementation of emergency plans. Furthermore, considering that M r reduces the variability of the between-event residuals 41 with respect to M w , we foresee its possible application in the development of ground motion prediction equations for seismic hazard assessment.

Method
Calibration of the empirical relationships IV2 S -E r and PD S -M 0 . Following previous studies 19 , we have derived an empirical relationship between IV2 S and E r . To this purpose, we considered 229 earthquakes, and we computed the theoretical seismic radiated energy following a procedure 52 based on the spectral integration of the squared theoretical Brune's velocity spectrum for S-waves 25 . These have been derived using the corner frequency (f c ) and seismic moment (M 0 ) computed in previous studies 8,53 .
Hence, we related the IV2 S estimated for each recording to E r radiated by the source assuming a linear model, where the attenuation of average IV2 S amplitude as a function of distance is expressed without assuming any a-priori functional form (i.e., non-parametric, or data-driven, approach): where the hypocentral distance R H range is discretized into N bin ; the index j = 1,…, N bin indicates the j-th node selected such that R H is between the distances r j ≤ R H < r j+1 ; the attenuation function is linearized between nodes r j and r j+1 using the weights w, computed as w j = (r j+1 − R H )/(r j+1 − r j ). The R H range 5-100 km is discretized into 19 bins with equal width (i.e., 5 km) (Fig. 5a). The minimum distance is fixed to 5 km given the lack of recordings at shorter distances. The coefficients A, B, C j are determined by solving the over-determined linear system (8) in a least-square sense. To fix the trade-off between A and C j , the attenuation is constrained to zero at r 2 = 10 km (i.e., at the upper boundary of the first distance bin). However, it is worth noting that changing the node constrained to zero corresponds to changes of both the A and C j coefficients which compensate each other, so that the coefficient B in Eq. (8) is unaltered.
Similarly, the relationship between the peak displacement PD S and M 0 has been expressed in the form: S H j j j j 0 1 The coefficients A, B, C j and D, F, G j are reported in Table S1. The calculated regressions have a R 2 equal to 0.84 and 0.87, respectively.
The scaling with distance of the derived attenuation models is shown in Fig. 5a and b, where the distance distribution of the observed IV2 S and PD S values corrected for the source scaling A + B log(E r ) and D + F log(M 0 ) are compared with the attenuation models C j and G j , respectively (i.e., where the index j indicates different bins of hypocentral distances). For both cases, the good agreement with the corrected data confirms the suitability of the obtained attenuation models. Figure 5c and d compare the observed IV2 S and PD S values corrected for attenuation effects, as modelled through the C j and G j coefficients, with E r and M 0 . The source scaling (black line) defined by the coefficients A, B and D, F capture well the trend in the data over the entire energy and moment ranges, as indicated by the agreement with the median value obtained for each event by averaging over more than 10 recording (yellow circles in Fig. 5c and d). Finally, Fig. 5e and f summarize the residual distributions between predicted and observed log(E r ) and log(M 0 ). Considering the values computed for each recording (gray), the residuals are unbiased with standard deviations equals to 0.56 and 0.29 for E r and M 0 , respectively; when the average values per event are considered (yellow), the standard deviations reduce to 0.27 and 0.1, respectively. Data and Resources. The waveforms used in this study have been obtained from European Integrated Data