Source Parameters of Moderate-To-Large Chinese Earthquakes From the Time Evolution of P-Wave Peak Displacement on Strong Motion Recordings

In this work we propose and apply a straightforward methodology for the automatic characterization of the extended earthquake source, based on the progressive measurement of the P-wave displacement amplitude at the available stations deployed around the source. Specifically, we averaged the P-wave peak displacement measurements among all the available stations and corrected the observed amplitude for distance attenuation effect to build the logarithm of amplitude vs. time function, named LPDT curve. The curves have an exponential growth shape, with an initial increase and a final plateau level. By analyzing and modelling the LPDT curves, the information about earthquake rupture process and earthquake magnitude can be obtained. We applied this method to the Chinese strong motion data from 2007 to 2015 with Ms ranging between 4 and 8. We used a refined model to reproduce the shape of the curves and different source models based on magnitude to infer the source-related parameters for the study dataset. Our study shows that the plateau level of LPDT curves has a clear scaling with magnitude, with no saturation effect for large events. By assuming a rupture velocity of 0.9 Vs, we found a consistent self-similar, constant stress drop scaling law for earthquakes in China with stress drop mainly distributed at a lower level (0.2 MPa) and a higher level (3.7 MPa). The derived relation between the magnitude and rupture length may be feasible for real-time applications of Earthquake Early Warning systems.


INTRODUCTION
The characterization of the seismic source in terms of earthquake magnitude and source radius (or length of the rupture) is now a routinely operation in any standard seismological laboratory. However, both parameters are generally computed off-line, through fairly complex procedures, mainly performed in the frequency domain. The seismic moment, for example, is estimated from the low frequency amplitude of displacement spectra. The source radius is typically obtained from the spectral corner frequency (Brune, 1970;Madariaga, 1976) or from time-domain, source duration measurements, generally available several minutes after the earthquake occurrence (Boatwright, 1980;Duputel et al., 2012). Although the fitting of spectral shapes is a straightforward operation, a major issue is the adequate correction of the observed spectra for path attenuation and site response effects.
In the context of Earthquake Early Warning (EEW), the point source characterization is usually based on the measurement of a few parameters (typically peak amplitude and/or characteristic period) in the early portion of the recorded P waves (3-4 s). These parameters are related to the earthquake size or to the peak ground shaking through empirical relationships (Allen and Kanamori, 2003;Kanamori, 2005;Wu and Zhao, 2006;Zollo et al., 2006;Böse et al., 2007;Wu and Kanamori, 2008;Zollo et al., 2010).
More recently, new strategies have been proposed to improve the accuracy of source parameter estimation for EEW applications and provide an estimate of the rupture area extent and the slip distribution on the fault. Among these strategies, some of them are based on the rapid inversion of geodetic and/or accelerometer data or fitting the spectrum in realtime (Allen and Ziv, 2011;Ohta et al., 2011;Colombelli et al., 2013;Caprio et al., 2011;Ziv and Lior, 2016). However, the rapid inversion of geodetic and/or accelerometer data need a catalog of the active faults for the construction of a rupture model plane in real-time (Colombelli et al., 2013). The azimuth dependency and the simplifying assumptions (e.g., directivity and segmentation) may introduce large discrepancies between modeled and observed spectra, leading to large variability in corner frequency estimates during the real-time spectrum inversion method (Ziv and Lior, 2016).
A second class of algorithms is based on the real-time spatial assessment of ground motion values. The FinDer algorithm (Böse et al., 2012), indeed, provides an estimate of the fault rupture extent and strike by continuously monitoring the spatial distribution of ground motions in real time. The results of this approach can be affected by the accuracy of near/far source classification and by the shortcomings in the GMPE used for the template generation (Böse et al., 2012(Böse et al., , 2015. More innovative approaches to EEW bypass the real-time source parameter estimation (location and magnitude) and use physics-based data assimilation techniques to directly predict the incoming evolution of ground shaking (PLUM) algorithm (Hoshiba and Aoki, 2015).
Recently, Colombelli and Zollo (2015) looked at the time evolution of the early P-wave information and used it as a proxy for the rupture process of earthquakes to extract the seismic moment and rupture extent of moderate-to-large Japanese earthquake records., Nazeri et al. (2019) explored a similar approach using strong-motion data of the 2016-2017 Central Italy sequence and estimated moment magnitude, fault length and average stress drop for each single event. The proposed method accounts for the effects of azimuth and distances by averaging the distance-corrected peak displacement among many stations, distributed over azimuth and distance, to approximate the moment rate function (MRF) and can be applied to general geometries with no need of prior knowledge of fault information. The proposed method is a remarkably simple and straightforward approach that does not require any complex calculation for automatically estimating two main source parameters (the earthquake magnitude and the expected length of the rupture) before the rupture has finished.
Since 2005, the National Strong Motion Observation Network System of China was established. Stations in the network were mainly equipped with force balance accelerometers with a broadband frequency range of 0-80 Hz, which ensures that the network is reliable to record large quantities of high-quality strong-motion data for research purposes (Li et al., 2008).
Following the idea of Colombelli and Zollo (2015) and Nazeri et al. (2019), in this study we use a database of moderate-to-large Chinese earthquake records and explore a similar approach to estimate the earthquake magnitude and rupture length, and to provide an approximate estimate of the average stress drop to be used for Earthquake Early Warning and rapid response purposes. The aim is twofold: 1) establish source scaling relationships for moderate to large earthquakes in China and 2) build the foundation for further studying the feasibility of a networkbased EEW method based on the time evolution of the early P-wave peak displacement amplitude.

Data Selection and Construction of LPDT Curves
For the present analysis, we selected the earthquakes occurred in China in the period 2007-2015. The magnitude of all the events (surface wave magnitude, Ms) varies between 4.0 and 8.0. To avoid the inclusion of bad quality data in our analysis, we selected seismic records with an epicentral distance smaller than 120 km, but for the M8 event we expanded the limit to 200 km and required that each event had at least three records. A total of 1293, 3-component accelerometric waveforms, relative to 88 earthquakes and 540 stations were used for the regression of P-wave peak displacement amplitude attenuation relationship, among which we selected 31 earthquakes (for a total of 617 records) with at least ten recording stations for the computation of the LPDT curve. Figure 1A shows the epicentral position of the selected earthquakes and the location of stations, in which two main seismic regions (Sichuan-Yunnan and Xinjiang regions) in China have been enlarged for clarity. Figure 1B shows the histogram distribution of the analysed records as a function of the epicentral distance and magnitude.
We identified the onset of the P wave on the vertical component of acceleration records, using a standard shortterm/long-term average method for automatic picking (Allen, 1978). Then, we visually inspected all the available waveforms and made manual picks where necessary, to adjust potential mistakes from the automatic picking algorithm. After removing the mean value and the linear trend, the acceleration waveforms are integrated once to velocity and twice to get displacement. Finally, we applied a 0.075 Hz high-pass Butterworth filter to remove the low frequency drift on displacement records. We impose the zero-crossing of the signal amplitude at the onset of the P-wave, to eliminate any potential residual noise contaminations resulting from the double integration operation. We then measure the absolute maximum of the initial P-wave amplitude on the vertical component of displacement (named Pd) using an expanding time window, starting at the arrival of the P-wave and moving forward with a time step of 0.01 s. The peak amplitude is related to the earthquake magnitude (M) and to the source-to-receiver distance (R) through an attenuation Frontiers in Earth Science | www.frontiersin.org March 2021 | Volume 9 | Article 616229 3 relationship of the general form (Wu and Zhao, 2006;Zollo et al., 2006): where Pd is the P-wave peak measurements and A, B and C are coefficients empirically determined from select dataset using 2 or 3 s after the P-wave arrival. Unlike the previous study, here we calibrated the coefficients by performing a specific least-squares multiple regression analysis, in which we fixed the distance attenuation coefficient (C) and chose a fixed length of the P-wave time window for the parameter measurements, which was set at 3 s for M ≤ 7 and 9 s for 5.5 < M ≤ 7, respectively. Further details about the estimation of the coefficients of Eq. 1 are provided in Supplementary Material. For the computation of the LPDT curve of each event, the peak amplitudes Pd of all records are measured at every P-wave time window and the distance-corrected amplitudes (logP c d) are obtained as logP c d logPd − ClogR. In order to avoid the contamination by the S-waves on the selected portion of the P-wave, we picked out four stations with clear seismic phase randomly from the dataset in each distance bin (every 20 km). A total number of 40 records were used, and their S-wave arrival time were manually picked to estimate the coefficients of the following equation: where Tp is the P-wave onset time, Ts is the arrival time of the S-wave, R is the hypocentral distance in km, b 0.13 is the coefficient derived from a linear regression analysis. To minimize any potential S-wave contamination, we regarded Tp+0.8bR as the expected S-wave arrival. As the time window increases, the stations with the expected S-wave arrivals were automatically excluded to make sure only P-wave part involved in the computation. Finally, the LPDT curve is obtained by averaging the distance-corrected amplitude of all the valid stations at each time window. The computation of the curves stops when the number of stations is less than a minimum of data (Five stations).
Observation and Modelling of LPDT Curve Figure 2 shows an example of the generated LPDT curve for the M4.6 event. The LPDT curve has an exponential growth shape with an initial increase, a gradual intermediate curvature and a final plateau level. Generally, the LPDT curve of larger event needs more time to reach the plateau and the plateau level of the curves scale with the final magnitude ( Figure 3A). The shape of the LPDT curve, as obtained from the average of many stations distributed over azimuth and distance, can be interpreted as a proxy of the Moment Rate Function (MRF), from the initial time up to its maximum peak value. Therefore, two essential features of the MRF, i.e., peak value and peak time, which are both related to the source properties, should be embodied in the LPDT curves. Following the idea of Colombelli and Zollo (2015), for near-triangular source time functions, the peak value of the MRF (related to the magnitude) will correspond to the plateau level of the LPDT curve, and the peak time of the MRF (related to rupture half-duration) is a proxy for the time at which the LPDT curve reaches its plateau level (Plateau Time). With this in mind, the magnitude and rupture duration can be estimated from the plateau level and the plateau time of the curves.
To model the LPDT curves, we fit data using the following function (Colombelli et al., 2020): where y 0 is fixed as the first point of the curve, P L is the interval between y 0 and the plateau level, a is the weighting factor which is set to 0.5, T 1 , and T 2 are the time parameters (here we define the larger value as T 2 ). T 1 controls the very initial part which usually has a faster increasing speed and T 2 represents the second part, whose increasing speed gradually becomes slower. This double corner time, exponential model accounts finely for the two different behaviors of LPDT curve-that is to say, a sharp increase to the plateau (ramp-like) for small events and a more gentle and smooth increase (exponential) for large events. The model parameters are shown in Figure 2.
We used a non-linear, weighted-fitting approach to model our curves, accounting for the standard error on each point of LPDT curves. Specifically, at each time step, the weight is obtained as: where SE is the standard error in each P-wave time window, N is the number of stations used for that time window. Figure 3B shows that the LPDT curves of all the events are quite well FIGURE 2 | The LPDT curve of M4.6 event and the model parameters.
The grey triangles stand for the corrected amplitude Pd C for the available stations at each P-wave time window, then the circle represents the average Pd C value of all these stations. Our exponential fit model is shown in black line. The time parameters T 1 and T 2 are shown in the empty triangles.
Frontiers in Earth Science | www.frontiersin.org March 2021 | Volume 9 | Article 616229 4 reproduced by the fitting model, with an average residual of 0.3. Generally, at the beginning of the curve computation, several stations from a broad range of distances and azimuths are involved in the calculation, so that the scatter of data is large and the fitting procedure gives a smaller weight to this part, as compared to the plateau of the curves, which is instead, well reproduced. A slightly larger (about 0.7) difference between the real data and the model is observed for the initial part of the curve of the M 8.0 Wenchuan earthquake, which could be related to the complexity of the source process of this peculiar event. Indeed, when looking at the seismic moment release of this event ( Figure 4), a small peak value is observed at 4-5 s, before the

RESULTS
We fit the LPDT curves with our exponential model and obtain the three relevant parameters mentioned above (T 1 , T 2 , P L ) while y 0 has been fixed to the first point amplitude. For simplicity, we defined a new variable called P L * P L + y 0 to represent the true plateau level of the LPDT curves. Both the amplitude parameter (P L *) and the two characteristic times (T 1 and T 2 ) scale with earthquake magnitude. Figure 5 shows the plateau level P L * as a function of magnitude. A good correlation between P L * and magnitude (the correlation coefficient reaches 0.95) can be found. The parameters T 1 and T 2 extracted from our fitting model are shown in Figure 6 as a function of magnitude. As the best fitting line indicates, both parameters linearly increase with magnitude (in logarithmic scale). Due to the very rapid initial increase of the curves, T 1 and T 2 are close to each other for small events, while they gradually separate when the magnitude becomes larger and the initial part of the curves increases gently.

Magnitude Estimation
Once the observed amplitude has been corrected for the distance effect, the LPDT value at each P-wave time window can be associated to a corresponding magnitude using the coefficients of Eq. 1 (Supplementary Table S1). Figure 3 shows the dynamic process of estimating magnitude for the LPDT curve. The y axis on the left stands for the distance corrected Pd, and the corresponding magnitude or estimated magnitude scale is shown on the right. The magnitude then can be estimated accurately based on the true plateau level of the LPDT curves and the coefficients of Eq. 1 listed in Supplementary Table S1. As shown in Figure 3, the occurrence of the plateau for large events (M > 6-7) needs more than 9 s after the P-wave arrival, suggesting that the typical approaches for the magnitude estimate using fixed 3-4 s P-wave time windows (PTWs) would provide underestimated magnitudes for such large events. Moreover, the B coefficient (calibrated using a 3 s PTW) could not be suitable to compute the corresponding magnitude based on the P L * obtained in a longer time window. We therefore choose two magnitude ranges, with two different PTW lengths, to calibrate and use the optimal coefficients (Supplementary Table S1) for magnitude estimation. In this way, when an event reached its plateau within 4 s, we use the relation coefficients A and B for fixed 3 s PTW, while we used the coefficients A and B established with a fixed 9 s PTW when its LPDT curve keeps increasing after 4 s. The estimated final magnitude based on the P L * for the LPDT curve and obtained with the two sets of coefficients for small and large events respectively, is plotted in Figure 7 as a function of the catalog magnitude. As it can be seen in Figure 7, most of the points are distributed around the dashed line representing the 1:1 relationship between the estimated magnitude and the catalog magnitude. The scatter of data is rather small, with an average estimated error of 0.23 magnitude units.

Prediction of Rupture Length and Estimation of Stress Drop Δσ
The rupture duration is the total duration of a seismic event, given by the whole-time length of the moment rate function (Vallée,FIGURE 6 | Scaling relationship between T 1 and T 2 as a function of magnitude. Triangles and circles represent T 2 and T 1 for LPDT curves, respectively. The error bars which represent the standard deviation of the parameters fitted by our exponential model are shown as black lines (where visible). The best fitting line for each parameter is shown as dashed line.
FIGURE 7 | The estimation of the magnitude and catalog magnitude. The dashed line corresponds to the 1:1 scaling. The triangles and circles indicate the estimated magnitude with P L * using the relation coefficients A and B for a fixed 3 s PTW and a fixed 9 s PTW, respectively.
Frontiers in Earth Science | www.frontiersin.org March 2021 | Volume 9 | Article 616229 6 2013). It is generally observed that the rupture duration scales with magnitude and is related to the rupture length, assuming a value for the rupture velocity (Wells and Coppersmith, 1994). As the moment rate function can be simplified as a symmetric triangle (Bilek et al., 2004), the peak time of the MRF, corresponding to the plateau time (T PL ) of the LPDT curve, is a measure of the Half-duration of each event (Colombelli and Zollo, 2015;Nazeri et al., 2019). Having in mind the Sato & Hirasawa model (Sato and Hirasawa, 1973) here we investigate the relation between parameter T 2 of LPDT curves and T PL , the time at which occurs the peak of the MRF. In the Sato & Hirasawa model, the rupture spreads radially outwards at a constant velocity with a circular fault, and stress drop (Δσ) and rupture velocity (Vr) are two relevant parameters controlling the earthquake rupture process. Since Δσ and Vr of earthquakes can vary significantly for each event (Allmann and Shearer, 2009), we performed a set of dedicated simulations to explore stress drop values between 0.05 and 20 MPa with rupture velocities between 0.5 Vs (S-wave velocity) and 0.9 Vs, for a total of 55 combinations of the two parameters.
For each given magnitude, we fixe Δσ and Vr, and generate the corresponding Sato & Hirasawa moment rate Function (SHF) by changing the polar angle of the observation point from 0 o to 90 o and computing the average SHF (an example of M 5 event shown in the Figure 8A). We then compute the log of the SHF and keep unchanged after the SHF reaches its peak to get a curve (hereinafter LSHF) with a similar shape of our LPDT curve. Since most of the selected earthquakes in our database occurred at an average depth of about 10 km, we set the Vp 6.2 km/s, Vs 3.4 km/s based on the velocity model for this region . Examples of the average SHF with fixed Vr 0.7 Vs and Δσ 0.1 MPa is shown in Figure 8, while examples with other Δσ values are shown in the Supplementary Material. For each available couple of stress drop and rupture velocity, we used the exponential model (Eq. 3) to fit the LSHF curve and extract the T 2 parameter for different magnitudes. As expected, we found that T 2 has linear relationship with T PL obtained directly from the peak time of the generated SHF (in logarithmic scale) for the entire magnitude range with a small deviation when exploring the Δσ and Vr, suggesting that the T 2 parameter extracted from the observed curves can be used to predict T PL : log 10 (T PL ) 1.111(0.051)log 10 (T 2 ) + 0.542(0.030) (5) For the circular model with a symmetric triangular shaped MRF, the obtained T PL can be regard as the Half-Duration (HD) of the source function. Considering all azimuthal coverage around the fault, the relation between averaged half-duration Frontiers in Earth Science | www.frontiersin.org March 2021 | Volume 9 | Article 616229 and the source radius was given as follows (Aki and Richards, 2002): where a is the source radius, Vr is rupture velocity, θ is the polar angle of the observation point and Vp is the P-wave velocity. Given the half-duration of the source, the source radius of the analysed events can be estimated. Figure 9 shows the predicted source radius as a function of magnitude and its corresponding half-duration with a fixed rupture velocity of 0.9 Vs. Based on the computation of stress drop (Δσ) for the circular model using the source radius (a) and the seismic moment (M 0 ) (Keilis-Borok, 1959), Δσ 7 16 M 0 a 3 , the theoretical scaling lines of the source radius as a function of M with a constant Δσ 0.1 and 10 MPa were given in the same figure as a comparison. The predicted source radius shows a similar increasing trend with the theoretical lines, indicating that the source radius of the analysed events has a consistent selfsimilar, constant stress drop scaling with magnitude. We fit the source radius with a weight-based fitting approach (same as Eq. 4, here, the SE is the standard error of the source radius computed from the predicted T PL and its error obtained by the error propagation theory) to obtain the best-fit constant stress drop of 0.4 MPa. In addition, we repeated the process by setting Vr 0.7 Vs and Vr 0.8 Vs, and found that the mean value of Δσ are 1.0 and 0.6 MPa, respectively.
The definition of rupture length changes in case of a circular fault (rupture length rupture radius) and long-rectangular fault (rupture length larger rectangular fault dimension). For a circular fault, the rupture lengths (source radius) predicted from LPDT curve were shown in Figure 10. As the figure indicates, most of our predicted rupture lengths as a function of magnitude agrees with the magnitude-rupture scaling relation studied by Cheng et al. (2019) Wells and Coppersmith (1994). However, our predicted results were smaller than the rupture length of two models in the moderate-large magnitude range (M 7.0 Lushan earthquake and M 8.0 Wenchuan earthquake). To better assess the performances of our LPDT scheme in predicting rupture length, we collect the rupture lengths of these two largest events derived from other methodologies. Our approach provides a rupture length of 34 km for the M 7.0 Lushan earthquake, which is consistent with the results of relocated aftershocks by  . The circular model is well suitable for small-to-moderate magnitudes (i.e., M < 6). Here, we also investigated the results of large magnitudes when applying the rectangular source model of Haskell. The far-field displacement radiated by a Haskell type fault model is equivalent to the convolution of two box-car functions of different amplitude and durations: rise time (τ) and rupture time (T R ). The resulting function has a trapezoidal shape with total duration given by the sum of τ + T R . The rupture time T R depends on the finite length of the fault (L) and the azimuth (θ) between source and receiver (Haskell, 1964): The rise time τ is independent of azimuth (Hwang et al., 2011) and can be obtained using the following relationship derived by Melgar and Hayes (2017) from a database of finite faults: log 10 (τ) −5.323 + 0.293log 10 (M 0 ) (8) Figure 11 shows an example of the average total duration of Haskell model with azimuth (θ) changed from 0 o to 180 o . Assuming that T PL is the middle point of the average trapezoid plateau, the following relationship between L and T PL can be obtained and be used for estimating rupture length of large events when assuming the Haskell model: We computed the predicted rupture length for the two largest events using the Haskell model. As shown in Figure 10, with the same T PL , applying the Haskell model appears to predict larger rupture lengths than the results using circular model. Specifically, we obtain an estimate of rupture length of 326 km without underestimation for Wenchuan earthquake, whereas the rupture length of Lushan earthquake has been overestimated (L 77.3 km). The circular model produced a triangle-like MRF, as the Sato & Hirasawa moment rate Function of M 5 event shown in the Figure 8A indicates, the peak time (T PL ) is almost half-duration. However, as the MRF of M 8.0 Wenchuan earthquake (download from USGS, shown in Figure 4) indicates, the MRF of some large events usually have the major peak occurring at the beginning of the rupture and followed by a long-time duration coda (Vallée, 2013). In this case, the circular model with a triangle-shape MRF may be not able to correctly reproduce this kind of source time function, and considering the peak time as the half-duration of the large events may result into underestimated rupture lengths. The averaged MRF from Haskell model shown in Figure 11 reaches its peak (T PL ) around 1/5 of the total duration, which is more like the characteristic of MRF of the large magnitude earthquakes stated above. Moreover, for large earthquakes, the width (W) of the fault is already saturated, i.e., equal to the thickness of the brittle fracturing zone in the lithosphere. Accordingly, the growth of the fault area with increasing seismic moment (M 0 ) is in the length direction only Cheng et al., 2019). Therefore, in the situation of L >> W for large events, the Haskell model is more appropriate. Since the rupture surface of Lushan earthquake presents the focal characteristics of smallmoderate earthquakes (L ≈ W) (Chen et al., 2013;Hao et al., 2013), the Haskell model is not suitable for this event leading to an overestimated rupture length.
According to the derived rupture length of the analysed events, the stress drop can be computed based on different geometrical faults. Madariaga (1977) proposed that the general form Δσ where W is the rupture width. Referring to the observation of Cheng et al. (2019) that large events with M > 6.7 have a nearly constant rupture width of ∼20 km, here, we set W to 20 km.  Li et al. (2017). As Figure 12B shows, the stress drops vary from 0.01 to 10 MPa. We divide the events into two groups based on Δσ 0.6 MPa and calculate the averaged stress drop for these two groups. The higher stress drop group (Δσ > 0.6 MPa) has an averaged stress drop of 3.7 MPa, which is close to the world-wide measured median value of 4 MPa (Allmann and Shearer, 2009). The lower stress drop group with an averaged stress drop of 0.2 MPa consists of 19 events, in which half of the events are the aftershocks of the M 8.0 Wenchuan earthquake and the M 7.0 Lushan earthquake. Consistent with Wang et al. (2018), the lower stress drop of aftershocks may result from the remaining locked parts on the fault plane of the mainshock. All these source parameters of the analyzed events (M ≥ 5.5) are summarized in Table 1.

DISCUSSION AND CONCLUSION
In this study, we generalized the approach proposed by Colombelli and Zollo (2015) to estimate the source parameters of a set of Chinese earthquakes with magnitude ranging from 4 to 8. The methodology is based on the use of the time evolution of the P wave (LPDT curve) as a proxy for the source time function to extract earthquake magnitude and rupture duration.
Comparing with the previous works by Colombelli and Zollo (2015) and Nazeri et al. (2019), we used the double corner time, exponential model proposed by Colombelli et al. (2020) for better modeling the behavior of LPDT curves. We improved the magnitude estimation based on the plateau level of LPDT curve, by using two different scaling coefficients with fixed C, which have been properly calibrated. Based on the analysis of 31 events in the magnitude range between 4 and 8, we found that the plateau level of LPDT curves has a strong correlation with magnitude (the coefficient of correlation is up to 0.95). Comparing with the catalog magnitude of the analyzed events, our predicted magnitude from the displacement data shows an average deviation of 0.23 magnitude units. The time at which these estimates are available implicitly depends on the event itself. Small events have a rapid initial increase and reach the plateau quickly, and for these events, the magnitude estimates would be available within a short time from the first P-wave detection. For the large events, like Wenchuan earthquake (M 8.0), the time needed to reach the plateau is longer. The near-source stations with a short fraction portion of the P-wave time-series cannot provide robust magnitude estimates for these large events, which is consistent with our observation that the LPDT curve at this stage is still growing and the plateau has not yet been reached. However, the expansion of hypocentral distance allows to capture longer P-wave signal portion of large events, and thus we provide a magnitude estimation of 8.3 for the Wenchuan earthquake after the LPDT curve reaches its plateau. By capturing the curve that reaches its plateau or looking at larger distances, our approach can estimate the magnitude of large events accurately, without saturation effects, typically observed when using shorter P-wave time windows (Lomax and Michelini, 2009;Bormann and Saul, 2009;Colombelli et al., 2012).
Together with the plateau level, the plateau time (T PL ) of the LPDT curve has also a clear scaling with magnitude, being approximately related to the half-duration of the source. In order to estimate the plateau time, we performed a series of simulations based on the Sato and Hirasawa (1973) circular rupture model and on the assumption that T PL corresponds to the peak time of the MRF. We generated a set of MRFs exploring Δσ and Vr values and studied the relation between T PL and the time characteristic parameter T 2 . Finally, we established a linear scaling relationship between the two parameters for predicting the T PL . Considering the circular model with a symmetric triangular-shaped MRF, the obtained T PL can be regarded as the Half-Duration (HD) of the events to predict the source radius. The obtained source radius in this study shows a consistent self-similar, constant stress drop scaling with magnitude. We obtained the best-fit stress drop (0.4 MPa) for the 31 analysed events, with fixing rupture velocity to 0.9 Vs. This value is lower than world-wide measured median value of 4 MPa (Allmann and Shearer, 2009), but it is comparable with the mean value of 0.5 MPa by studying the strong-motion recordings of the Wenchuan aftershocks (2008-2013) (Wang et al., 2018).
One major result of this paper concerns the determination of the scaling law of earthquakes in China. We obtained the rupture length of the analysed events using circular model and found the M 8.0 Wenchuan earthquake has a shorter predicted rupture length with comparing to other results. We realize that for the largest events in the sequence, the circular fault model may underestimate the total rupture length. Thus, we suggested to estimate the rupture length of the large events (M > 7) assuming the Haskell, rectangular fault model. Our predicted rupture length of different source models as a function of M is close to the rupture scaling relationship proposed by Cheng et al. (2019) and Wells and Coppersmith (1994).
The most critical issue of the proposed methodology is to assume a triangular shape to represent the source function. The assumption of a circular fault model can be inadequate to describe the complexity of the source for large events, for which multiple peaks of moment release can occur. Moreover, the rupture process of an earthquake is a chaotic process, and the simple exponential LPDT curve may hardly represent the complex nature of this process, especially at the very beginning. However, the now massive and well documented observation of LPDT empirical curves analysed in many worldwide seismic regions and in a wide magnitude and distance range confirm that the general ramp-like behavior of the curves is a universal evidence (Colombelli et al., 2012;Colombelli et al., 2015;Colombelli et al., 2020;Melgar et al., 2015;Trugman et al., 2019) which is well matched by the chosen exponential model (Colombelli et al., 2020). When the curve has reached its peak, the schematic representation with two parameters (plateau time and amplitude), may properly catch the essential features of the process we are interested in.
Concerning the anelastic attenuation effect on LPDT curves, the overall features of the curves are essentially dominated by the event magnitude, which mainly controls the plateau level and the time to reach this plateau, for each curve. Using similar LPDT curves, Colombelli et al., (2020) recently investigated the effect of anelastic attenuation on the shape of the curves. This was done both theoretically and through a careful and detailed analysis of real earthquake data. They found that the chosen anelastic attenuation parametrization (Eq. 1) is appropriate to correct the distance attenuation effect on LPDT curves and this does not produce a systematic effect on plateau and T PL parameter determination in all the observed curves in the same distance range.
The source parameters of the larger events are a major focus of -and motivation for -the proposed methodology. Although the selected dataset has a relatively wide magnitude range from 4.0 to 8.0, there is an obvious lack of large-magnitude data for events with magnitudes between 7 and 8. The proposed methodology would be better verified through the application to other seismic areas, once that the proper calibration of scaling coefficients has been done. It is worth to mention that the findings (e.g., rupture scaling relation and stress drop distribution) obtained through applying the methodology to the selected events (31 earthquakes) are mainly based on the data from Sichuan-Yunnan region, and may not be applicable to other areas as they are.
We notice that in the presence of a jump in the increasing shape of the LPDT curve, our approach likely fits the curve with a longer T 2 thus resulting in a lower stress drop value. Supplementary Figure S2 shows the LPDT curve for the M 6.6 Ludian earthquake: a jump in the curve is clearly visible around 3.5 s. A similar magnitude event (M 6.4) which has a normal LPDT curve shape without jumps is selected for comparison. Our double corner time, exponential model (Eq. 3) fits the LPDT curve of this M 6.4 event with a T 2 0.66 s. However, due to the presence of jump, our exponential model tends to find a parameter (T 1 and T 2 ) with slower growth to fit the middle part (2 ∼ 4 s) of the LPDT curve for M 6.6 Ludian earthquake. This resulted in our model fitting the LPDT curve with a longer T 2 1.8 s. We further justify the effect of LPDT curve shape on the exponential fitting model in Supplementary Figure S3. The virtual LPDT curve with normal LPDT curve shape was generated as follow: the initial increasing part (0 ∼ 1.85 s) and the approach plateau part (4.26 ∼ 6 s) are the same as the observed LPDT curve, and only the middle part, we add points with exponential growth. After using our exponential model to fit the curves separately, we obtained a T 2 1.46 s for the virtual LPDT curve with a normal shape and a slightly longer T 2 (1.8 s) for the observed LPDT curve with a jump in the middle. Thus, based on the relationship between T 2 and T PL , the predicted T PL (plateau time) of LPDT curve with normal shape then will be 5.3 s, whereas the T PL of the LPDT curve with jump will be 6.7 s. The jumps on the curve can be given by different effects, both related to the source process itself, such as multiple peaks on the Moment Rate Function, and/or to propagation effects, such as the arrival of intermediate phases (preceding the S-wave arrival), arriving at a receiver within the coda of P-wave phase. However, the presence of jumps on the LPDT curves is not frequently observed in our catalogue, only happens in two events. Meanwhile, the jump in the LPDT curve does not have a serious impact on the derived parameters. Thus, we believe that, on average, this is a negligible effect on the overall estimation of rupture length and stress drop.
The shape of the curves may change in real-time, when we do not have all the available stations. As a perspective of the work, we could evaluate the feasibility of application in real-time, that can be relevant for EEW applications. The accurate estimation of the rupture extent and magnitude at the early stage of the process can be a useful piece of information to add to the early shake-map computation, and the estimation of ground shaking can be strongly improved by considering the effect of the finite fault (Colombelli et al., 2013). While in this off-line study we used the post-earthquake location instead of the real-time estimation, a reliable estimation of the earthquake location is needed for the real-time application. Having in mind that the far-field stations must wait enough time to reach the plateau of the curve for large events, we need more time to get the plateau information. Hence, the real-time application performance and the timeliness based on the network distribution for this approach will be further studied. A possible method could be that we can estimate the final curve at each time step with a given probability and study the realtime curve reached how many percent of the final curve (maybe after T 2 ) can give a reliable probability for estimating final curve.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
YW applied the methodology and drafted the manuscript; SC and AZ critically revised the paper and supervised the study; JS and SL provided suggestions on the results. All authors contributed to the article and approved the submitted version.