Earthquake Magnitude Estimation from High-rate GNSS Data: A Case Study of the 2021 Mw 7.3 Maduo Earthquake

: Peak ground displacement (PGD) and peak ground velocity (PGV) are critical parameters during earthquake early warning, as they can provide rapid magnitude estimation before rupture end. In this study, we used the high-rate Global Navigation Satellite System (GNSS) data from 55 continuous stations to estimate the magnitude of the 2021 Maduo earthquake in western China. We used the relative positioning method and variometric approach to acquire real-time GNSS displacement and velocity waveforms, respectively. The results showed the amplitude of displacement and velocity waveforms gradually decreased with increasing hypocentral distance. Our results showed that the ﬂuctuation of PGD magnitudes over time is smaller than that of PGV magnitudes. Nonetheless, the earthquake magnitudes estimated from both methods were consistent with their counterparts (Mw 7.3) reported by the United States Geological Survey (USGS). The ﬁnal magnitude estimated from the PGD and PGV methods were Mw 7.25 and Mw 7.31, respectively. In addition, our results highlighted how the number of high-rate GNSS stations could inﬂuence the stability and convergence time of magnitude estimation. conﬁrmed applicability reliability real-time


Introduction
According to the China Earthquake Network Center, at 18:04:13 (UTC) on 21 May 2021, an Mw 7.3 earthquake occurred in Maduo County, Qinghai Province, China. The epicenter was located at (98.34 • E, 34.59 • N) and had a focal depth of 17 km (Figure 1). Field investigations immediately after the event [1,2], and Interferometric Synthetic Aperture Radar (InSAR) observation obtained within a few days after the main shock [3][4][5] confirmed that the seismogenic fault of the Maduo earthquake was the Kunlunshankou-Jiangcuo Fault, a secondary fault~70-80 km to the south of the East Kunlun Fault within the Bayan Har block. This unexpected earthquake was the largest earthquake to have occurred in China since the 2008 Mw 7.9 Wenchuan earthquake [6][7][8], challenging the conventional perspective that the Bayan Har block acts as a quasi-stable block.
Over the past 20 years, high-rate (≥1 Hz) Global Navigation Satellite System (GNSS) stations have been set up across mainland China to study crustal deformation. During the Maduo earthquake, the densely distributed high-rate GNSS stations (Figure 1) around the epicenter recorded ground displacement and velocity waveforms, providing valuable data to investigate the source process of the earthquake. Moreover, these data can help in testing algorithms, such as the rapid estimation of magnitude, in the geodetic-based earthquake early warning (EEW) system that is under development [9]. Real-time processing of GNSS data is challenging but vital for EEW. Currently, there are three categories of methods for real-time GNSS kinematic positioning: Precise Point Positioning (PPP) [10], relative positioning [11] and the variometric approach [12]. PPP relies on precise orbit and clock products, and requires a relatively long convergence or re-convergence time (20-30 min) due to complicated error corrections and parameter estimation [13,14]. Relative positioning has an accuracy of sub-centimeters, but it needs at least one stable reference station to acquire displacements and thereby might be affected by the reference station [15,16]. Nonetheless, it has been widely used to extract real-time coseismic displacements. The variometric approach uses broadcast ephemeris to obtain GNSS velocities with an accuracy of mm/s in real time [12]. This method does not need reference stations and could shorten the convergence time when compared with PPP. Site velocities could be easily integrated to displacements with an accuracy of cm, but drifts may occur during the process, resulting in error accumulation [17,18].
Magnitude is a key parameter in EEW. In estimating earthquake magnitude, traditional seismic data, such as ground acceleration recorded by strong-motion accelerometers and ground velocity recorded by broadband seismometers, may suffer from some limitations. For instance, broadband seismometers are prone to clip and go off-scale when recording in the near-field of large earthquakes (Mw > 7.0) [19]. This leads to magnitude saturation, a well-documented condition in large earthquakes whereby an EEW system underestimates the true event magnitude. Although acceleration records have the advantage of amplitude unsaturation in the near-field of large earthquakes, one needs to integrate acceleration to displacement or velocity to calculate the magnitude. Due to factors such as instrument tilt and rotation, the integration would suffer from baseline drift during this process [20].
GNSS is capable of measuring long periods down to the static offset (0 Hz) and could, therefore, ameliorate some of the above limitations in determining earthquake magnitude. In the time domain, although seismometric methods only take a few seconds (0.5-4 s) for the P wave to estimate the magnitude of an event, their predictions typically saturate for M 7+ earthquakes [21]. GNSS could account for static surface displacement accumulated with the arrival of the S wave, enabling the estimation of finite-fault slip and the event's moment magnitude [22].
To investigate the applicability of GNSS in magnitude estimation of large earthquakes in the Tibetan Plateau, in this study, we focused on the recent Maduo earthquake. When using GNSS to calculate earthquake magnitude, peak ground displacement (PGD) is often used. The regression model between PGD and earthquake magnitude was first established using the seismogeodetic observations from five earthquakes (Mw 5.3-9.0) in Japan and California [23]. The regression model was tested and updated using more earthquake cases in recent years [24][25][26]. The past few years benefited from the accumulation of high-rate GNSS data in earthquake case studies, and a regression model between earthquake magnitude and peak ground velocity (PGV) of high-rate GNSS was established [16]. Using the data from 22 earthquakes with magnitudes ranging from Mw 6.0 to Mw 9.1, Fang et al. [16] found that the PGV magnitudes were comparable with PGD magnitudes.
In this study, we collected high-rate GNSS data from the 2021 Mw 7.3 Maduo earthquake. First, we used the relative positioning method to obtain the coseismic displacement waveforms and the variometric approach to obtain the velocity time series at high-rate GNSS stations. PGD and PGV magnitudes were then estimated for each station. Finally, we compared the differences between PGD and PGV magnitudes and discussed their implications in geodetic-based EEW.

GNSS Data and Data Processing
High-rate GNSS stations within 550 km of the epicenter recorded the coseismic displacement and velocity waveforms of the 2021 Mw 7.3 Maduo earthquake. We collected the high-rate (1 Hz) GNSS data from 55 continuous stations (Figure 1), 14 of which were from the Crustal Movement Observation Network of China (CMONOC). The remaining were from the Qinghai Continuously Operating Reference Stations (QHCORS) network.
To acquire displacement waveform results at each high-rate GNSS station, we used the relative positioning method implemented as the TRACK module in the GAMIT/GLOBK software [27]. During processing, we used the International GNSS Service (IGS) final orbits products and absolute antenna phase center model. We selected the SCYY station, which is located in the Sichuan Province and is about 850 km away from the epicenter, as the reference station. To avoid masked coseismic signals, we did not perform the postprocessing of spatial filtering of displacement waveforms. To confirm that our processing results were not affected by the reference station, we used the PPP method [28] to process the high-rate data of the reference station, and the result ( Figure S1) showed that the reference station was not affected by coseismic deformation.
To obtain site velocity time series results, we used the variometric approach implemented as the SNIVEL software package [18] to process the high-rate GNSS data. This approach uses real-time broadcast ephemeris data to estimate the high-precision real-time velocities at GNSS stations [12,18]; it uses dual-frequency data to form a linear combination of L1 and L2 narrow lanes [18]. The above practice reduces noise level and corrects ionospheric error and tropospheric delay. No external data, such as IGS products, were needed during data processing.

PGD and PGV Magnitudes
The three-component displacement time series was used to extract PGD value from each GNSS station, which is given as: where N d (t), E d (t) and U d (t) represent the displacement components of the north, east and up in meters, respectively. The regression model between PGD and moment magnitude is given as [23][24][25][26]: where Mw is the moment magnitude; R is the hypocentral distance in km; and A, B and C are the regression coefficients. We used the regression coefficients (A = −5.919, B = 1.009 and C = −0.145) from Ruhl et al. [26], which were estimated from GNSS observational data taken during 29 moderate-to-large (Mw 6.0-9.0) earthquakes. At each high-rate GNSS station, we calculated the moment magnitude at each epoch until the moment magnitude reached convergence. We then averaged the moment magnitudes from all stations to obtain the final moment magnitude of the Maduo earthquake.
Using the three-component velocity time series, the PGV value at each GNSS station was calculated as: where N v (t), E v (t) and U v (t) show the velocity components of the north, east and up in m/s, respectively. The regression model between PGV and moment magnitude is expressed as [16]: where Mw is the moment magnitude; R is the hypocentral distance in km; and A, B and C are the regression coefficients. We used the regression coefficients (A = −5.025, B = 0.741 and C = −0.111) from Fang et al. [16], which were extracted from modeling of GNSS observational data taken during 22 moderate-to-large (Mw 6.0-9.1) earthquakes. Similar to the PGD magnitude, the PGV magnitude was updated at every epoch until it reached convergence. The average magnitude at all GNSS stations was taken as the final moment magnitude of the Maduo earthquake.  Figure 1. Figure S2 shows the displacement waveforms at all high-rate GNSS stations. We notice that vertical displacements were less sensitive to Maduo earthquake than horizontal displacements. This is to be expected, given the strike-slip focal mechanism in the 2021 Maduo earthquake. Therefore, we focus on horizontal displacements in the following section. The amplitude of displacement waveforms gradually decreases with increasing hypocentral distance. Moreover, as the hypocentral distance increases, permanent horizontal coseismic displacement decreases.   Figure S3 shows the velocity waveforms at all high-rate GNSS stations. It is obvious that before the arrival of the S wave, velocity waveforms at all stations were nearly flat. After the earthquake, all stations experienced obvious velocity fluctuations. The amplitude of velocity waveforms at GNSS stations attenuates as the hypocentral distance increases. Similar to displacements, vertical velocities were less sensitive to seismic waveform than horizontal velocities. Nonetheless, quite a few GNSS stations observed significant vertical velocity fluctuations.

PGD/PGV Results and Magnitudes
We extracted PGD and PGV values at each high-rate GNSS station from the displacement and velocity waveforms, respectively. Figure 3 shows PGD and PGV values as a function of hypocentral distance. It is obvious that PGD and PGV values decrease with hypocentral distance increasing. PGD and PGV values fluctuate near the line, showing the theoretical relationship between PGD/PGV and hypocentral distance in Mw 7.3 earthquakes. Figure 4 shows the PGD magnitude evolution. We obtained a maximum and minimum value of PGD magnitude of Mw 8.0 and Mw 6.8, respectively. These results indicate that, at different high-rate GNSS stations, PGD magnitude varies to a sizable extent, which might be related to site effects and/or radiation pattern. More reliable magnitude estimates can be obtained by averaging the results at multiple GNSS stations [29]. The final average magnitude was estimated at Mw 7.25, which is slightly smaller than the moment magnitude (Mw 7.3) reported by USGS, indicating that the PGD magnitude from GNSS is reasonable.  Slightly different from the PGD magnitudes, PGV magnitudes show larger deviation. The maximum PGV magnitude was estimated at Mw 8.7, whereas the minimum was Mw 6.6. The average magnitude after convergence was Mw 7.31, slightly larger than that from PGD estimation but consistent with the magnitude (Mw 7.3) reported by USGS.
The above results show that, in spite of different magnitude estimates at different GNSS stations, both PGD and PGV can provide plausible earthquake magnitude estimates that are consistent with the magnitude derived from the seismological method. We present more analysis in the following section regarding the time efficiency of PGD and PGV methods in estimating magnitude.

Discussion
PGD and PGV are critical parameters during EEW; they can provide rapid magnitude estimation before rupture end [30][31][32], thereby issuing reliable warning information to the public. During the Maduo earthquake, because the majority of high-rate GNSS stations were far from the epicenter and unevenly distributed, we discuss how the number of high-rate GNSS stations might influence the stability and convergence time of magnitude estimation. In addition, we provide outlooks for the application of PGD and PGV magnitude estimation methods. Figure 5a shows the PGD and PGV magnitudes as function of the number of GNSS stations. The number of GNSS stations increases as the epicenter distance increases. Overall, PGD and PGV magnitudes fluctuate around Mw 7.3, with maximum deviations of Mw 0.21 and Mw 0.22, respectively. Compared with the PGV magnitudes, PGD magnitudes are more stable and fluctuate less. However, the final PGV magnitude (Mw 7.31) is closer to Mw 7.3 than the PGD magnitude (Mw 7.25). In particular, we notice that in the case of four stations, both PGD and PGV underestimated the magnitude. When the number of GNSS stations increases to about eight, the PGD method can obtain a stable magnitude, while the PGV method requires more stations (>15) to obtain a stable magnitude. In this sense, the PGD method seems to be more suitable for applications in areas where GNSS stations are sparse. We chose to use the eight stations within 170 km of hypocentral distance to discuss the sensitivity of magnitude convergence time to the number of high-rate GNSS stations. Figure 5b,c show the convergence time of PGD and PGV magnitudes as a function of the number of GNSS stations. There are two main features. First, in the case of the same number of stations, the convergence time for PGD and PGV magnitudes was almost identical. For example, when the number of GNSS stations was four, the convergence times for PGD and PGV magnitudes were all at 50 s after the earthquake onset, but the first-alert magnitude by the PGD method (Mw 5.5) was significantly higher than that of the PGV method (Mw 3.1). Second, with the increases in the number of GNSS stations, the time for the magnitude to converge also increases. For instance, when the number of GNSS stations was eight, the convergence times for PGD and PGV magnitudes were at 73 s after the earthquake onset; this is to be expected, because different stations received S waves at different times.
In terms of time efficiency, we discuss the PGD magnitude shown in Figure 5b. Using four stations from the current GNSS array, we obtained the first alert 6 s after the earthquake onset; 50 s later, the magnitude reached convergence with estimate at Mw 7.1. When the number of GNSS stations was eight, the convergent magnitude (Mw 7.5) was at 73 s after the earthquake onset. The more stations there were, the more accurate the final estimated magnitude would be. We note that these convergence times are longer than the~40 s of the source rupture. This is because almost all the GNSS stations were in the far-field of the Maduo earthquake. Nonetheless, the GNSS displacement time series data provide irreplaceable constraints for rapid finite-fault slip inversion.
EEW is currently being carried out on a global scale [9,21,33]. PGD and PGV from high-rate GNSS data provide a way to quickly estimate the earthquake magnitude [16,32], as shown during the Maduo earthquake in this study. Although this method does not require complex fault rupture models, there are some issues that need attention and resolution. Firstly, compared with the relative positioning method that requires at least one reference station, the variometric approach uses broadcast ephemeris to obtain GNSS site velocities and is thus more suitable for early warning networks where stations are relatively sparse. However, as we emphasized in the Introduction, integrating velocity to displacement may cause error accumulation. Therefore, it is necessary to develop methods of removing errors effectively. Secondly, compared with seismic data with a high sampling rate (≥50 Hz), the sampling rate of GNSS data (1-10 Hz) is significantly lower. Therefore, when using GNSS data to estimate the earthquake magnitude, it is inevitable that high-frequency information would be lost, which might result in deviations in earthquake magnitude estimate. Consequently, it is imperative to develop methods of effectively integrating geodetic data with seismic data [34][35][36]. Finally, compared with PGD magnitude, PGV magnitude at different stations fluctuates greatly, meaning that the regression model between PGV and moment magnitude may still need to be improved when dealing with earthquakes in Tibet. Regarding the results of our work, integrating PGD and PGV magnitudes seems to be a plausible method of magnitude estimation.

Conclusions
In this study, we collected the high-rate GNSS data from 55 stations around the 2021 Mw 7.3 Maduo earthquake epicenter. Using the relative positioning method to process these data, we obtained displacement time series at each GNSS station. In addition, we used the variometric approach to process GNSS data and obtained velocity time series at each GNSS station. Displacement and velocity data were used to extract PGD and PGV values, which were then used to estimate PGD and PGV magnitudes.
Our data processing results showed the amplitude of displacement and velocity waveforms gradually decreased with increasing hypocentral distance. PGD and PGV fluctuated near the line (Mw 7.3), which reflects the relationship between earthquake magnitude and hypocentral distance. PGD and PGV magnitudes were consistent with their counterparts (Mw 7.3) from USGS, with absolute deviation of 0.05 and 0.01 magnitude units.
Our results confirmed the applicability and reliability of real-time earthquake magnitude estimation from high-rate GNSS data.