Validation of Peak Ground Velocities Recorded on Very-high rate GNSS Against NGA-West2 Ground Motion Models

Observationsofstronggroundmotionduringlargeearthquakesaregenerallymadewithstrong-motion accelerometers. These observations have a critical role in early warning systems, seismic engineering, source physics studies, basin and site amplification, and macroseismic intensity estimation. In this manuscript, we present a new observation of strong ground motion made with very high rate (>= 5 Hz) Global Navigation Satellite System (GNSS) derived velocities. We demonstrate that velocity observations recorded on GNSS instruments are consistent with existing ground motion models and macroseismic intensity observations. We find that the ground motion predictions using existing NGA-West2 models match our observed peak ground velocities with a median log total residual of 0.03-0.33 and standard deviation of 0.72-0.79, and are statistically significant following normality testing. We finish by deriving a Ground Motion Model for peak ground velocity from GNSS and find a total residual standard deviation 0.58, which can be improved by ~2% when considering a simple correction for Vs30.


Introduction
Ground motion models (GMMs), traditionally known as ground motion prediction equations (GMPEs), are empirical relationships between an earthquake source and the ground motion expected at a station. GMMs will commonly incorporate magnitude, distance attenuation, site amplification terms, and source terms to determine the peak ground accelerations, velocities, displacements, and spectral accelerations/displacements at different periods (PGA, PGV, PGD, and SA/SD, respectively). The coefficients in GMMs are determined empirically using observations from real earthquakes on strong-motion accelerometers or broadband seismometers if they are unclipped. GMMs are utilized in engineering seismology for building code design (Bommer et al., 2010;Katsanos et al., 2010), in ShakeMap and PAGER (Prompt Assessment of Global Earthquakes for Response) generation for rapid assessment of impact and ground motions after large earthquakes (Wald et al., 1999;Allen et al., 2009), in paleoseismic studies (Rasanen et al., 2021), and in earthquake early warning systems (Allen and Melgar, 2019;Meier, 2017;Thakoor et al., 2019). Many GMMs are derived for specific regions where epistemic uncertainties may be consistent across the region, but some use generalized global datasets or only include one style of faulting. Each GMM will have a parametric specific range in which the model is valid, i.e. a distance or magnitude limit, and most generally do not include very large magnitude events or far-field observations. There is also a question of whether or not PGV observations for large ground motions recorded with inertial sensors are recording the true ground motions due to sensor rotations and tilts (Boore et al., 2002;Clinton, 2004). In this study, we provide a new observation which can faithfully record high PGV, Global Navigation Satellite System (GNSS) derived velocities.
For earthquake and tsunami early warning systems, the importance of GNSS displacements for large earthquake characteristization has been demonstrated (Crowell et al., 2009(Crowell et al., , 2012(Crowell et al., , 2016Blewitt et al., 2006;Grapenthin et al., 2014;Minson et al., 2014;Murray et al., 2018;Williamson et al., 2020). The addition of this data is critical to properly characterizing the impacts of large earthquakes (M > 7) where traditional seismic methods begin to saturate (Melgar et al., 2013a). GNSS displacements do not saturate since they are computed in a non-inertial reference frame, that of the fixed center of the Earth. The GNSS observations capture both the coseismic static offsets, important for understanding moment release and total slip, and dynamic motions, which can be leveraged for rapid magnitude determination and kinematic inversions. However, issues with phase ambiguity fixing, cycle slips, and loss of satellite lock can lead to large errors that obscure ground displacements when computing GNSS displacements in real-time. These processing issues can lead to displacement excursions of many meters, a major problem for real-time analysis in early warning systems. GNSS displacements are also noisier than seismic observations due to path errors between the satellites and receivers, most notably the tropospheric and ionospheric delays (Melgar et al., 2020). One potential solution is to compute velocities rather than displacements at the GNSS stations. This approach is commonly referred to as the 'variometric approach', and it consists of performing a single difference in time between satellite orbital positions and the raw GNSS phase observables (Colosimo et al., 2011;Benedetti et al., 2014;Geng et al., 2016;Grapenthin et al., 2018;Shu et al., 2018;Crowell, 2021). The variometric approach is more sensitive to smaller magnitude ground motions and when integrated into displacement, leads to lower noise than traditional positioning since the various errors in real-time positioning do not change appreciably over short time periods (Shu et al., 2018;Dittmann et al., 2022a); in this study, we will show that GNSS velocities can be resolved for earthquakes as low as M4.9 in the near-field (< 25 km), over a full magnitude unit less than what is possible with real-time GNSS displacements (Melgar et al., 2020;Goldberg et al., 2021).
In this manuscript, we first describe how the full GNSS velocity dataset, spanning a magnitude range from 4.9-9.1, is generated and provide statistics of the data. Next, we use 3 NGA-West2 GMMs and published finite fault models where available to predict ground motions for all station-event pairs and compare the predictions against our peak ground velocity observations. We also look at the statistics for specific earthquakes and conditions, and test for normality in the residuals using a Lilliefors test. Finally, we generate a GMM with magnitude and distance scaling terms and solve for the coefficients using two different distance measures.

Data and Processing Methods
The variometric approach for geophysical applications was first presented by Colosimo et al. (2011) for geo-physical applications. In this method, a single time difference is performed on the GNSS phase observables and on the orbital positions. For a single frequency, L 1 , and satellite, s, the simplified variometric observation equation is at a receiver, r, (1) ∆L 1,s,r = ∆ρ s,r + c(∆τ r − ∆τ s ) + ∆T s,r − ∆I 1,s,r + ∆M s,r + ∆B 1,s,r + ∆e 1,s,r where ρ s,r is the range between satellite and receiver, c is the speed of light, τ r and τ s are the receiver and satellite clock biases respectively, T s,r is the tropospheric component between satellite and receiver, I 1,s,r is the ionospheric delay on the L 1 frequency between satellite and receiver, M 1,s,r is the L 1 multipath component for a given satellite and receiver, B 1,s,r is the ambiguity term on the L 1 frequency for a given satellite and receiver, and e 1,s,r are any uncharacterized noise sources at the L 1 frequency. The ∆ indicates that we are taking a difference between the observations at the current time, t, and the observations at the prior time step, t−1. Within Equation 1, many of the terms drop out as they typically do not appreciably change over small time periods, notably the multipath, the satellite clock errors, the troposphere, the ionosphere, and the fractional cycle and integer ambiguities, assuming no cycle slips are present. Receiver clock drifts do need to be solved since the drift rate can exceed the nanosecond per second level, which would equate to tens of centimeters per second.
To solve for receiver velocities, we use the SNIVEL (Satellite Navigation derived Instantaneous Velocities; Crowell (2021)) software package, which currently only uses the Global Positioning System (GPS). SNIVEL forms the narrow-lane (N L) combination of L 1 and L 2 , which is defined as where f L1 and f L2 are the frequencies of the L 1 and L 2 bands, 1575.42 MHz and 1227.60 MHz, respectively. The narrow-lane combination has an effective wavelength of 10.7 cm, which is less than the wavelengths of L 1 (19.0 cm) or L 2 (24.4 cm) alone, which slightly reduces the noise on the derived velocities. For determining the orbital velocities, we use the GPS broadcast orbits, which are well modeled and readily available in real-time. To solve for receiver velocity, the basic observation equation in Equation 1 is Taylor expanded about the range, ρ, which allows for the generation of a design matrix based upon the direction cosines between satellite and receiver. We then set up the following linear observation model at a given receiver for n satellites: (3) where x s,r , y s,r , and z s,r are the Cartesian coordinates of the satellite and receiver, and v is the velocity of the receiver. We make two simple phase corrections to account for subtle atmospheric signal path variations using the readily available in real-time hydrostatic tropospheric correction of Niell (Niell, 1996) and the Klobuchar ionospheric correction (Klobuchar, 1987). The velocities and clock drift rate can be solved through ordinary least squares. We additionally weight the least squares problem through a diagonal matrix of elevation angle between satellite and receiver since lower elevation observations will have more noise.

SNIVEL Dataset
The data that we used was recorded by six networks (number of earthquakes followed by observations in parentheses), the Plate Boundary Observatory (27; 178), COCOnet (10; 46), TLALOCnet (5; 27), Geonet New Zealand (9; 158), the TU-CWU (Tribhuvan University -Central Washington University) network in Nepal (4; 29), and the Italian RING network (6; 71). A map view of the distribution of events is shown in Figure 1. The first three networks were historically separate, however, they now comprise the federated Network of the Americas (NOTA) and their data is archived at UNAVCO (Murray et al., 2019). Since 2007, the standard practice at UNAVCO has been to download any available 5 Hz GNSS data within a given radius after a significant event within North and Central America. The earliest NOTA-recorded earthquake we include in this study is the 2009 Mw 7.3 Honduras earthquake. Geonet, operated by the Institute of Geological and Nuclear Sciences (GNS), has downloaded 10 Hz raw GNSS data after significant events since 2013. The 2016 Mw 7.8 Kaikoura earthquake is the best recorded event in our dataset, with 122 observations and PGV values greater than 40 cm/s at 8 sites and greater than 95 cm/s at 2 sites. The TU-CWU network was able to record at 5 Hz the 2015 Mw 7.8 Gorhka earthquake in Nepal, the Mw 7.3 aftershock two weeks later and two additional aftershocks. Three of the sites during the Gorhka earthquake exceeded 50 cm/s. The RING network operated by Istituto Nazionale di Geofisica e Vulcanologia (INGV) consolidated data from several smaller networks within Italy (ISPRA, DPC, Regione Lazio, Regione Abruzzo, Leica ITALPOS, and Topcon NETGEO) and recorded most stations at 10 Hz and several at 20 Hz for events during the Amatrice-Norcia sequence in 2016 and the Emilia-Romagna sequence in 2012. A 20 Hz recording at station ARQT during the October 30, 2016 Mw 6.6 Amatrice earthquake exceeded 100 cm/s and 5 additional stations exceeded 30 cm/s. In total, we have 509 observations from 61 earthquakes ranging in magnitude from 4.9 to 8.2. Figure 2 shows the distribution of these observations as a function of rupture distance and an overview of all the events with the PGV values at every station is provided (see Data Availability). We low pass filter all of the waveforms to one-quarter of the sampling rate, 1.25 Hz for the 5 Hz observations, 2.5 Hz for the 10 Hz observations, and 5 Hz for the 20 Hz observations using a 4-pole zero-phase Butterworth filter. All waveforms are visually inspected to ensure the peak velocity is related to shaking rather than noise. We do find that for the lower magnitude events (M < 6), the low pass filtering of the waveforms does reduce the peak velocity values, which is predicted when looking at the corner frequency dependence on magnitude (e.g., Joyner, 1984). When using an average stress drop of 30 bars and a shear wave velocity of 3 km/s, the corner frequency of a M5 earthquake would be roughly 0.6 Hz. This value drops precipitously, and by M8, the corner frequency is roughly 0.02 Hz. For consistency, we retain the low pass filtering as we find it beneficial to reduce the high frequency noise within the GNSS velocities as they are plagued by higher order wet tropospheric, ionospheric, and multipath noise. The peak value of each component is taken separately and the station PGV value is the maximum of the three individual components (north, east, or vertical); note that this is the same approach taken by the GMMs. For smaller levels of shaking, the vertical component was unable to record any ground motion and these observations are excluded from our dataset (the vertical component is on average 3-5 times noisier than the horizontal components). Figure 2 shows the PGV values versus mean rupture distance for the entire dataset.

Seismogeodetic Dataset
In addition to the SNIVEL dataset, we include observations from two earthquakes that were processed through a multi-rate Kalman filter, the 2011 M9.1 Tohoku-oki and 2003 M8.2 Tokachi-oki earthquakes in Japan on the Geonet network operated by GSI (Crowell et al., 2009;Bock et al., 2011;Melgar et al., 2013b). For both of these earthquakes, 1-Hz GNSS displacements were computed and then Kalman filtered with collocated strong-motion accelerometer data from either K-NET or Kik-net. Since this dataset includes a combination of both seismic and geodetic data, it is termed 'seismogeodetic' herein. The Kalman filter effectively integrates the accelerometer under the displacement constraint of the GNSS and produces both velocity and displacement waveforms at the sampling rate of the accelerometer. Melgar et al. (2013a) showed that this method was better able to reconcile ground motions in the low frequency end than simply integrating accelerometer data. In total, 174 observations for the two Japanese events at 100 Hz are included. Similarly to the SNIVEL dataset, we low-pass filter the velocity observations to 25 Hz.

Ground Motion Models Used
We use three of the NGA-West2 (Next Generation Attenuation for Western United States, 2.0) ground motion models (GMM) to compare our observed PGV values: Chiou and Youngs (2014), Boore et al. (2014), and Campbell and Bozorgnia (2014), herein referred to as CY14, BSSA14, and CB14 respectively. The NGA-West2 database consists of events between magnitudes of 3 and 7.9, and most of its observations are between 0 and 400 km . The distance distribution between our dataset and NGA-West2 is comparable, Figure 1 Map view of all earthquakes in this study. The solid black circles are the epicentral locations of the earthquakes and the red circles are sized by magnitude.
although we are more skewed towards the larger magnitude end since velocities cannot be recorded on GNSS instruments well below M5.5 ( Figure 2). We include observations out to rupture distances of~800 km.
For events that we do not have slip models for, we treat them as point sources, otherwise we use the complete distance descriptions within the GMMs. The slip models we use are the National Earthquake Information Center (NEIC) official finite fault models (e.g., Hayes, 2017) which we extract directly from the geoJSON files provided on the USGS event overview pages (e.g., https://earthquake.usgs.gov/earthquakes/eventpage/ ci38457511/finite-fault). Between the different GMMs, there are 5 different distance measurements: Joyner-Boore distance (R jb ), rupture distance (R rup ), distance to the surface projection of the updip edge (R x ), depth to the updip edge of the rupture (Z tor ), and hypocentral depth. When treating the earthquake as a point source, R rup is the hypocentral distance, R jb and R x are the epicentral distance, and Z tor is the hypocentral depth. We also use the mean rupture distance (R p ) proposed by Thompson and Baltay (2018), which replaces the traditional rupture distance, R rup , with a slip scaled distance to account for the heterogeneity in slip distributions. For R p , we use a value of p = −2.0 for the power law weighting, which is the optimal value for PGV (Thompson and Baltay, 2018). While we make this correction, we do note that none of the GMMs were validated using mean rupture distance and we are simply using it as an additional comparison. We use the USGS global VS30 database to approximate the shear wave velocity in the upper 30 meters (Heath et al., 2020). Also, for basin terms (Z 1.0 and Z 2.5 ), we use the equations within the GMMs that relate VS30 to the basin terms. We disregard any hanging wall, directivity, and dip dependent terms within all of the GMMs for simplicity.

Results
The log residual between the GMM prediction and the observed PGV values, ln(P GV GM M ) − ln(P GV GN SS ), are shown as a function of mean rupture distance, R p , for the three GMMs in Figure 3. Table 1 shows the median and standard deviations for a number of scenarios using both R p and R rup (note that BSSA14 has the same performance for R rup and R p since it only uses R jb ). For comparison, the self-reported total residual standard deviation for the three NGA-West2 GMMs are 0.65, 0.54, and 0.58 for BSSA14, CY14, and CB14 respectively. For all three GMMs, the histograms of the log residuals are tightly distributed about zero. There are however differences between the three GMMs, mainly in the performance for the Kaikoura earthquake (grey shaded histograms on Figure 3). All three GMMs underestimate the observed PGV values for the Kaikoura earthquake, with BSSA14 performing the worst with a median residual of -1.29, and CY14 (median residual -0.45) performing slightly better than CB14 (median residual -0.73) when using R rup . Part of the underesti-

Figure 2
Distribution of GNSS velocity data used in this study. The upper left panel shows the PGV values versus the mean rupture distance, in log-log space. The mean rupture distance is defined in Thompson and Baltay (2018), and we use p = −2 here. The different symbols signify the sampling rate of the data. The histogram on the right shows the frequency of PGV observations. The histogram on the bottom shows the frequency of the mean rupture distance. mate for the Kaikoura earthquake is due to the strong northward directivity and the complexity of the rupture onto more than a dozen crustal faults (Hamling et al., 2017). The log residual using R rup for Kaikoura is better than R p primarily due to the shorter rupture distances making up for the impacts of directivity and shallow slip. When the Kaikoura earthquake is removed from the dataset, the log residual standard deviation dropped considerably, down to 0.61-0.67 when using R p . For BSSA14, the performance appreciably changes when removing Kaikoura, with it performing better that any of the other GMMs using R rup , which is unsurprising given the wide distribution of residuals shown in Figure 3b. No other event in the dataset has such a strong negative residual between the recorded ground motions and the GMMs. Surprisingly, all three GMMs model the Tohokuoki earthquake well, with most stations showing a residual less than 1 log unit even though none of the GMMs used ground motion data for events greater than Mw 8. We postulate that the rather compact nature of the Tohoku-oki rupture coupled with the further distances to stations allows for the GMMs to predict PGV well.
For all three of the GMMs, there is a slight underestimation of ground motions at further distances (> 400 km). This is unsurprising as this distance is outside the specified distance range for the GMMs. When we exclude stations further than 400 km, there is an improvement for BSSA14 and CB14, but no change for Figure 3 The log-residual between PGV from the GMMs and GNSS. Panels a, c, e, and g show the residuals as a function of mean rupture distance for GMMs BSSA14, CY14, CB14, and CDDG22 (this study) respectively. Panels b, d, f, and h show the histograms of the residuals for the GMMs to their left for all observations. On the histograms, the smaller grey bars in the foreground show the distribution for the 2016 Mw 7.8 Kaikoura earthquake.
CY14, indicating the CY14 is more sensitive to further rupture distances than the self-reported model limit of R rup < 300 km. We also computed statistics for two additional event subsets: subduction zone events and non -subduction events (minus Kaikoura). None of the three GMMs were not directly developed for the subduc-tion environment, so it is important to understand any systematic biases that may arise due to the tectonic environment. When discounting Kaikoura, there are 268 observations of subduction earthquakes and 293 observations from primarily strike-slip faults. When using R p , the median residuals for the non-subduction events were all negative, indicating the GMMs are underestimating ground motions. For the subduction events, the standard deviations were all lower than the nonsubduction events, and the median residuals were better for BSSA14 and CB14. This result is somewhat paradoxical in that the GMMs we are comparing against in this study were developed primarily for upper-crustal faults, so a better fit for those events would be expected, however, the shallower events will have far more variability with regards to source terms, distance measurements, and directivity such that a greater variability in the ground motion residuals would be observed. For subduction events, the source distances are generally greater so much of the GMM complexity can be averaged out. We would expect even better fits by using subduction zone specific GMMs (e.g., Parker et al. (2022)).

Normality Testing
While the distribution of results and statistics shown in Figure 3 and Table 1 are indicative of reasonable model and observation agreement, we test the null hypothesis that the log residuals of the GMMs are drawn from a normal distribution. In order to do this, we perform a Lilliefors test, which is similar to a Kolmogorov-Smirnov test except the data is allowed to be from any generalized normal distribution and the test statistic is more stringent (Lilliefors, 1967). The null hypothesis for the test is the log residuals are drawn from a normal distribution. For this test we choose a significance level of 0.05, which the p-value needs to be larger than, and the test statistic needs to be smaller than the critical value for the number of data points in the test. While there is no requirement that our data residuals be normally distributed, it does lead to higher confidence that any variations between our data and the GMMs are due to random Gaussian noise and not due to biases in the observations or GMMs. For 683 observation points, the Lilliefors test statistic critical value is 0.035. Of the GMMs tested, only CY14 passes the Lilliefors test and cannot reject the null hypothesis, with a p-value of 0.14 and a test statistic of 0.032. This is further illustrated when looking at the Quantile-Quantile (Q-Q) plots in Figure  4, where the theoretical quantiles, those drawn from a normal distribution, are plotted against the data/sample quantiles. While near the centers of the theoretical quantiles, both CB14 and BSSA14 perform well, at the edges there are considerable outliers. The Q-Q plot for CY14 is good out to 2.5 quantiles. From these results, we can confidently say that the GNSS velocities can be represented well by CY14 at a random normal level.

Noise Characteristics of GNSS Velocities
To characterize the relative noise levels on the GNSS velocities, we processed 30 minutes of data at station ARQT in the Italian RING network at 1, 5, 10, and 20 Hz on day 300 of 2016, between the M5.5 and M6.1 Norcia earthquakes. Table 2 shows the standard deviations for the three components of motion at the 4 sampling rates. Along all three components of motion, the standard deviation increases with higher sample rates, which has been shown in other studies (e.g., Shu et al., 2018), and the vertical component is roughly twice as large as the horizontal components; also, there is a strong autocorrelation of noise between the components and a regional correlation of noise due to the similar constellation geometries. The reason behind this is the errors in the orbits, multipath, and clocks are fairly constant at high frequency, but we are taking incrementally smaller and smaller time steps, which to first order, leads to an increase in noise. When looking at the power spectral content at these sample rates (Figure 5), we see that the noise behavior is more complex. We see between periods of 0.3 and 2 s that the power is progressively smaller for higher sample rates, however, there is an increase in power for smaller periods until the noise becomes roughly white below 0.25 s. Since for higher sampling rates, the data is in the higher frequency, higher noise part of the spectrum, the standard deviations are higher, but the higher sampling rates reduce the noise at given periods due to reducing the higher order noise in the GNSS observation model. This would indicate that it is beneficial to sample GNSS velocities at the highest possible sample rate, but then re-sample the data to 5-10 Hz, or wherever temporal aliasing will be minimized based upon the frequency content of the earthquake (e.g., Joyner, 1984). Indeed, Table 2 shows that when we re-sample the 20 Hz phase observations to 5 Hz and process through SNIVEL, we obtain a standard deviation of roughly half that obtained by processing the data directly at 5 Hz.

Ground Motion Model Development
We developed a preliminary GMM (herein referred to as CDDG22) using the same formalism as in Thompson and Baltay (2018) and Goldberg et al. (2021) where we only consider magnitude scaling (F M ) and distancescaling (F R ) terms such that The distance-scaling term is separated into an anelastic attenuation term (c R2 ) and a geometrical spreading/magnitude-dependent attenuation term The magnitude-scaling term roughly follows the form of Chiou and Youngs (2014), however we use the expression directly from Goldberg et al. (2021) that is simplified to ignore small magnitude terms The coefficients in Equations 4-6 are solved through a least-squares regression and can be found in Table 3. We test two different distance measures for R: R p and R rup . Figure 3 shows the log-residual plots of our GMM using R p and the observed GNSS PGV values. The standard deviation is 0.58 and 0.59 using R p and R rup respectively, and the residuals are much more tightly distributed about zero. We also do not see any appreciably anomalous earthquakes like Kaikoura when using the  Table 2 Standard deviation for GNSS velocities at station ARQT processed at 4 different sampling rates. The time period covered is between 17:30 and 18:00 UTC on October 26, 2016. σ has units of cm/s.

NGA-West2
GMMs (e.g., Figure 3). In CDDG22, we have not considered any site terms, such as Vs30. In Figure  6, we plotted the residual error for CDDG22 against Vs30 and see that there is a small trend that is linearly modeled as When we apply this site correction to our residuals, we reduce the standard deviations to 0.57 and 0.58 for R p and R rup respectively. While this is not an appreciable reduction (∼2%), it does indicate that further investigation into site-specific corrections is warranted as we expand out our dataset in the future. Statistically, the residuals between our GMM and the data pass the Lilliefors test (p-value = 0.067, test statistic 0.034) and the Q-Q plot (Figure 4) shows excellent agreement out to 2.5 quantiles.

Discussion
The standard mode of GNSS displacement positioning today is to either process the raw observations at a central processing center that can accumulate corrections -0.258 -0.303 Table 3 The coefficients for the CDDG22 GMM using either R p or R rup . The coefficients are defined in Equations 4-7.
for satellite orbits, clocks, fractional cycle biases, and tropospheric corrections or to perform the computations onsite while transmitting the corrections to the station. Within PPP computations, it is imperative to properly compute the integer cycle ambiguities before precise displacements can be obtained, which requires segregating the fractional cycle biases and integer ambiguities (e.g., Geng et al., 2012). This is generally an iterative process and requires many minutes for the solution to converge to its stable and precise estimate. Even after properly correcting for all sources of error, realtime displacement observations generally have a noise level of 1-2 cm in the horizontal components and 5 cm in the vertical (Melgar et al., 2020) over short time periods and long period drifts due to constellation geometry and multipath are evident at periods greater than 30 seconds. The GNSS velocity approach does not require any external corrections and can be easily deployed onboard GNSS receivers since it uses the broadcast orbital information from the satellites. The accuracy of our GNSS velocity estimates in the horizontal is roughly 0.6, 1.2, and 2.3 cm/s at 5, 10 and 20 Hz respectively, and

Figure 4
Quantile-quantile plots of the log residual between the GMM determined PGV and the GNSS derived PGV using R p . The quantiles assume a normal distribution with the data corrected and normalized to a mean of 0 and standard deviation of 1. The red dashed lines indicate a 1-to-1 correspondence for perfectly normal data.
these values could be further improved by using International GNSS Service Ultra-Rapid orbits and clocks and temporal re-sampling from higher rates (e.g., Shu et al., 2020).
Given the accuracy of the GNSS velocity estimates with respect to the GMMs, their use in ShakeMap generation should be considered. ShakeMaps are a powerful post-earthquake evaluation tool that provides information on ground motion (peak ground acceleration, velocity, spectral accelerations) and shaking intensity (in the form of Modified Mercalli Intensity, MMI). These maps are generated from a combination of GMMs relating an earthquake source to ground motions, instrumental recordings of ground motions, and com-munity intensity reports (i.e., Did You Feel It). The work presented here has demonstrated that GNSS velocities can be added to the instrumental recordings within ShakeMaps. This is important because the station distribution of seismic and GNSS networks are vastly different due to the motivating factors behind the installations. Indeed, Grapenthin et al. (2018) demonstrated for the 2017 Mw 7.1 Iniskin earthquake the enhanced value of GNSS derived velocities in a region with sparse seismic network coverage. Seismic networks were primarily installed either near well known seismic sources or near major population centers. GNSS networks were installed primarily to aid the surveying community and are thus more evenly distributed geo-   Table 3.
graphically. GNSS networks have also been targeted at high strain accumulation regions. This is now somewhat changing with collocating of instruments and expansion of seismic networks, but the different station distributions allow for us to provide observations of PGV in new locations, further constraining the interpolation schemes used in ShakeMap. To show the efficacy of GNSS velocities within ShakeMaps, we computed several sets of ShakeMaps for the 2016 Mw 6.6 Norcia, which are shown in Figure 7. For this demonstration, we replaced the input instrumental data with the PGV values recorded at the GNSS stations ( Figure   7a). We see in the near-field, there are considerably more GNSS stations than seismic stations (Figure 7b), which leads to considerably higher intensities. We also computed the ShakeMap using both the seismic observations and the GNSS velocities, which is shown in Figure 7c. This model is closer to the original ShakeMap, however, in the near-field, there is upwards of 0.8 MMI difference between the two models ( Figure 7d). While there are many issues still to explore with the optimal ways to incorporate this data into ShakeMaps regarding the appropriate weighting schemes, this demonstration shows a promising potential improvement in near-real-time intensity characterization.
While we have not explicitly shown the utility of GNSS derived velocities to earthquake early warning, there are several ways in which we envision these observations improving geodetic early warning. Many studies have shown the value of GNSS displacements to rapid magnitude and slip determination for large earthquakes, but the aforementioned issues with real-time positioning require clever logic within operational early warning systems. For example, the G-FAST system that is operating on ShakeAlert development servers, requires a gross outlier filter, a magnitude and time dependent uncertainty scheme, a time dependent minimum displacement filter (on 4 stations), and a minimum seismic magnitude filter (Murray et al., 2021). These logic filters built onto the back-end of G-FAST have the effect of throttling messages from the system from all events except those that we have the highest confidence in. Incorporating GNSS velocity streams into G-FAST could remove or reduce the levels of these filters by lending more confidence to the displacement streams. For example, ambiguity resolution issues on the displacement streams do not appear on the velocity streams, so we can use this stream to flag parts of the time series that there is no expected shaking. Moreover, the ground velocity observations can be used to appropriately window displacement time series or to select only those stations that should have a significant displacement signal and exclude those that are effectively noise, thus improving the precision and accuracy of the geodetic source models within G-FAST. For example, Dittmann et al. (2022b) trained a random forest classifier to select only parts of the GNSS velocity time series that had earthquake related ground shaking and showed a true positive rate of roughly 90% for earthquakes greater than M5 and out to hypocentral distances greater than 1000 km for larger events. Finally, GNSS derived velocities can be used directly in early warning systems either in existing seismic algorithms that rely on velocity observations or through magnitude scaling. Fang et al. (2020) showed a simple PGV scaling relationship, with a similar form to peak ground displacement scaling used in early warning (Crowell et al., 2013), to determine earthquake magnitude with uncertainty of 0.26 magnitude units.

Conclusions
We have demonstrated that GNSS derived velocities, recorded at 5 Hz or greater, are capable of characterizing strong ground motions for moderate to large earthquakes without going off scale. These observations agree well with the three NGA-West2 GMMs that we compared against, with the best agreement with Chiou and Youngs (2014) using R rup and with Campbell and Bozorgnia (2014) using R p . Deriving our own simple GMM directly from our PGV values and published USGS finite fault models yields roughly a 20% reduction in the log residual down to 0.58-0.59; this value is in line with the total residual reported by the three NGA-West2 GMMs using the seismic database, between 0.54 and 0.65 log units. Most importantly, our dataset includes PGV records from many large earthquakes, with almost half the observations (333) coming from M > 7.5 earthquakes, and this study provides true unfiltered records of strong ground motion in a non-inertial reference frame.

Acknowledgements
This work is funded by the US Geological Survey External Grants Program, grant number G22AP00038 to University of Washington. Additional funding is provided by USGS Cooperative Agreement G21AC10529 to University of Washington. This material is based on services provided by the GAGE Facility, operated by UNAVCO, Inc., with support from the National Science Foundation, the National Aeronautics and Space Administration, and the U.S. Geological Survey under NSF Cooperative Agreement EAR-1724794. We acknowledge the New Zealand GeoNet project and its sponsors EQC, GNS Science, LINZ, NEMA and MBIE for providing data/images used in this study. We also acknowledge RING GPS data by Istituto Nazionale di Geofisica e Vulcanologia, which are licensed under a Creative Commons Attribution 4.0 International License and outlined in the following reference: INGV RING Working Group (2016), Rete integrata Nazionale GNSS, doi:10.13127/RING. We would like to thank three anonymous reviewers who helped improve the manuscript.