A Non-ergodic Spectral Acceleration Ground Motion Model for California Developed with Random Vibration Theory

A new approach for creating a non-ergodic $PSA$ ground-motion model (GMM) is presented which account for the magnitude dependence of the non-ergodic effects. In this approach, the average $PSA$ scaling is controlled by an ergodic $PSA$ GMM, and the non-ergodic effects are captured with non-ergodic $PSA$ factors, which are the adjustment that needs to be applied to an ergodic $PSA$ GMM to incorporate the non-ergodic effects. The non-ergodic $PSA$ factors are based on $EAS$ non-ergodic effects and are converted to $PSA$ through Random Vibration Theory (RVT). The advantage of this approach is that it better captures the non-ergodic source, path, and site effects through the small magnitude earthquakes. Due to the linear properties of Fourier Transform, the $EAS$ non-ergodic effects of the small events can be applied directly to the large magnitude events. This is not the case for $PSA$, as response spectrum is controlled by a range of frequencies, making $PSA$ non-ergodic effects depended on the spectral shape which is magnitude dependent. Two $PSA$ non-ergodic GMMs are derived using the ASK14 and CY14 GMMs as backbone models, respectively. The non-ergodic $EAS$ effects are estimated with the LAK21 GMM. The RVT calculations are performed with the V75 peak factor model, the $D_{a0.05-0.85}$ estimate of AS96 for the ground-motion duration, and BT15 oscillator-duration model. The California subset of the NGAWest2 database is used for both models. The total aleatory standard deviation of the two non-ergodic $PSA$ GMMs is approximately $30$ to $35\%$ smaller than the total aleatory standard deviation of the corresponding ergodic $PSA$ GMMs. This reduction has a significant impact on hazard calculations at large return periods. In remote areas, far from stations and past events, the reduction of aleatory variability is accompanied by an increase of epistemic uncertainty.


Introduction
Ground-motion models (GMMs) are used to estimate the distribution of a ground-motion intensity measure (IM ) for a given earthquake scenario. The most common IM is pseudo-spectral acceleration (P SA) as it is a good estimator of seismic loading for a wide range of structures. P SA is defined as the absolute maximum response of a single-degree-of-freedom oscillator (SDOF) to an input ground motion. SDOFs are defined by their natural period (T 0 ) or natural frequency (f 0 = 1/T 0 ) and damping (ζ); in GMMs, typically, T 0 ranges from 0.01 to 10 sec and, ζ is equal to 5%. The response of the oscillator depends on the frequency content and timing (compactness of energy) of the ground motion. From the entire frequency content of the ground motion, the response of the oscillator mainly depends on the amplitudes of the frequencies near and below f 0 . Therefore, at small T 0 (high f 0 ), the response of the oscillator depends on the entire frequency content of the ground motion (i.e. spectral shape) and not just a narrow frequency bin. This makes the coefficients of a P SA GMM at small T 0 magnitude dependent even for linear effects, as the shape of spectral acceleration response spectrum changes with magnitudes. The peak of a spectral acceleration response spectrum will be at 0.1 sec for a magnitude (M ) 3 event and at 0.3 sec for a M 7.5 event ( Figure 1); this means that at small magnitudes, the P GA scaling (e.g. V S30 coefficient) will be consistent with the scaling of T 0 = 0.1 sec, while at large magnitudes, the P GA scaling will be consistent with T = 0.3sec. This is also observed by Stafford et al. (2017), who showed that the linear site amplification factors are magnitude and distance dependent. A detailed discussion the differences between the scaling of F AS and P SA is given by Bora et al. (2016). Most P SA GMMs do not explicitly account for the magnitude dependence of the coefficients, such as the V S30 scaling or distance scaling; instead, they often use a limited range of magnitudes where the magnitude dependence of the coefficients is not pronounced. For instance, the data-set that was used in the development of the NGA West1 GMMs had a limited set of magnitudes that ranged from M 4.5 to M 8 Power et al. (2008). The approach of using a smaller range of magnitudes works when developing an ergodic GMM, as there is enough number of moderateto-large magnitude events globally to estimate the coefficients, but it can be problematic when developing a non-ergodic GMM.
For the NGA West2 GMMs, the data set was extended to down to M 3 with the objective of setting the reference ergodic model that could be used to evaluate regional differences in the site, path, and source terms based on small magnitude data. The NGA West2 GMMs modified the magnitude scaling to capture the average effect of the magnitude dependence of the coefficients, but this does not accurately model the magnitude dependence of the site and path effects.
GMMs fall into two main categories: ergodic GMM and non-ergodic GMM. Ergodic GMMs assume that the statistical properties of a ground motion IM do not change in space (Anderson and Brune, 1999), and therefore, earthquakes and recordings from all around the world can be merged into a single dataset to estimate the GMM coefficients. Models developed under this assumption tend to have stable median estimates but large aleatory variability. Some models developed with the ergodic approach are: the NGA West GMMs for California Abrahamson et al. (2008), and the Douglas et al. (2014) GMM for Europe. Non-ergodic GMMs recognize that source, path, and site effects are systematically different at different parts of the world and account for these differences in the model development. Non-ergodic GMMs have smaller aleatory variability than ergodic GMMs, but in areas with sparse data, where the systematic effects are unknown, the reduced aleatory variability is accompanied by an increase in the epistemic uncertainty of the values of the median ground motion. The use of non-ergodic GMMs in Probabilistic Seismic Hazard Analysis (PSHA) is very promising, as the reduction in aleatory variability can have a large impact on the seismic hazard at large return periods, improving the accuracy of the site-specific hazard. A more in-depth discussion of ergodic and non-ergodic GMM is provided in the accompanying paper Lavrentiadis et al. (ress).
The estimation of the non-ergodic terms requires a large set of regional data. To achieve that, the datasets used in the development of non-ergodic GMM need to have a wider range of magnitudes to include the more frequent small-to-moderate earthquakes. It is this expansion of the magnitude range that makes the magnitude dependence of the GMM coefficients a more significant issue in non-ergodic GMMs. One solution to this problem is, first, develop a non-ergodic GMM for an IM whose scaling does not suffer from the magnitude dependence, as P SA does, and then for a scenario of interest, calculate the non-ergodic P SA based on the non-ergodic IM estimate.
The effective amplitude spectrum (EAS), defined in Goulet et al. (2018), is one such IM : the EAS is a smoothed rotation-independent average power Fourier amplitude spectrum (F AS) of the two horizontal components of an acceleration time history. In EAS, the amplitude at each frequency is independent of the amplitudes of the adjacent frequencies making the coefficients of an EAS GMM magnitude independent. Random vibration theory (RVT) provides a framework to calculate P SA from EAS. It relies on extreme-value statistics to estimate the peak response of the oscillator directly in the Fourier domain; it does not require a phase-angle spectrum to first convert the ground motion in the time domain to compute the peak oscillator response. RVT has been used in the past to compute P SA based on F AS from seismological theory (Hanks and McGuire, 1981;Boore, 1983Boore, , 2003 Other studies, such as Boore and Joyner (1984), Liu andPezeshk (1999) Bora et al. (2015) and, Boore and Thompson (2012), focused on semi-empirical adjustments to the RVT framework to correct for the assumptions not satisfied by ground motions, mainly the fact that acceleration time histories are not stationary signals. More recently, Kottke et al. (ress) used RVT to develop an ergodic P SA GMM for the eastern US based on an ergodic EAS GMM for the same region.
In this study, we developed two non-ergodic P SA GMM. The average P SA scaling is determined by backbone ergodic P SA GMMs. The non-ergodic effects are defined in terms of non-ergodic P SA factors which are estimated by combining the Lavrentiadis et al. (ress) non-ergodic EAS GMM with RVT.

Ground-Motion Data
A subset of the NGAWest2 data-set (Ancheta et al., 2014) was used in this study. The selected subset contains the earthquake and stations that are located in California, western Nevada, and northern Mexico. Recordings that were flagged as questionable in Abrahamson et al. (2014) were removed from the regression subset Figure 2 shows the spatial distribution of earthquakes and stations. Most of the stations are located in Los Angeles, Bay Area, and San Diego metropolitan areas, whereas spatial density of the stations is lower in less populated areas, such as northerneastern California. The regression data-set contains 7520 records from 185 earthquakes recorded at 1410 stations. Figure 3 shows the magnitude-distance distribution of the data and the number of records per frequency. The magnitude of the earthquakes ranges from 3.1 to 7.3, and the distance of most records ranges from 10 to 200 km. The usable frequency range of the majority of EAS records spans from 0.4 and 20Hz. The minimum usable frequency of most P SA records is 0.5 Hz.

Model development 3.1 Random-Vibration Theory
RVT uses Parseval's theorem and extreme-value statistics (EV S) to estimate the P SA based on the frequency content (i.e. F AS) and duration of a ground motion. Parseval's theorem is used to calculate the root-mean-square of the oscillator's response (x rms ) to the input ground motion, and a peak factor (P F ), based on EV S, is used to estimate the absolute peak response of the oscillator, which is the definition of P SA, based on x rms . P F s assume that the ground motion is a stationary stochastic process, and that it can be described as a band-limited white Gaussian noise with zero mean. The first assumption means that the amplitudes of the ground motion are identically distributed, and the second assumption means that the phase angles of the ground motion are   (b) number of P SA and EAS recordings per frequency used in the regression analysis randomly distributed. Although, earthquake ground motions violate both assumptions, numerous studies have shown that RVT provides P SA estimates that are in agreement with observed ground motions (Hanks and McGuire, 1981;Boore, 1983Boore, , 2003

Oscillator Response
The response of an oscillator to a ground motion can be computed by convolving the ground motion with the impulse response (IR) of the oscillator. The IR is the response of an oscillator to a very brief acceleration pulse; that is a Dirac delta function. For an SDOF oscillator, the Fourier transform of the impulse response is: where, f 0 is the natural frequency of the oscillator, and ζ is the damping of the oscillator. As an example, Figure 4 shows the P SA impulse response, in time and Fourier domain, for an SDOF oscillator with f 0 = 2Hz and ζ = 5%. In the Fourier domain, the convolution is performed by multiplying the ground motion's F AS with IR; therefore, the response of an SDOF oscillator to a ground motion is: The x rms of the oscillator's response is defined as: where D rms is a measure of the duration which is defined in Section 3.1.4. Parseval's theorem states that the amount of energy in the time domain is equal to the amount of energy in the Fourier domain ( with m 0 being the zeroth moment of F AS. The k th moment of F AS is defined as:

Peak Factor
The peak factor relates the x rms with the maximum response of the oscillator (x max ), which is the definition of the P SA.
In general, P F s fall into two main categories: those based on the Cartwright (1956) peak factor, abbreviated as CLH56, and those that are based on the Vanmarcke (1975) peak factor, abbreviated as V75.
In the first group, the CLH56 peak factor assumed that the peaks of a time history occur independently according to a Poisson process. In a series of papers, Boore and colleagues (Boore, 1983;Boore and Joyner, 1984;Boore, 2003) developed peak factors (BJ83) based on a reformulated version of CLH56 and removed an integrable singularity. Davenport (1964) proposed the a peak factor model (D64) based on an asymptotic form that approximates CLH56 for long time histories.
The main difference between V75 (Vanmarcke, 1975(Vanmarcke, , 1976) and the P F s of the first group is that V75 dropped the Poisson process assumption. Because of this, V75 P F accounts for the time spend outside the threshold, which is important for a narrow-band process, and considers that the peaks could be clustered in time, which is important for a wide-band process. Der Kiureghian (1980) noted that the D64 peak factor overestimates the number of zero crossings, and developed a new P F model (DK80) by modifying D65 P F so that it is asymptotically consistent with V75. V75 and D80 are in general agreement, but they deviate in time histories with a small number of zero crossings.
The V75 P F is selected for the development of the non-ergodic P SA GMM. V75 is preferred over the group of P F that are based on CLH56 due to the simplified assumptions in CLH56, and the complete form of V75 is preferred over the asymptotic forms, as the former is more accurate for the wide range of ground motions considered in this project. This choice is consistent with the P F used in Kottke et al. (ress).
V75 expressed the probability distribution of the peaks as a first-passage problem. For a Gaussian process, the first-passage probability (i.e. the probability of no crossing) a ±a threshold (type-D barrier) in the time interval (0, t) is equal to: where r is the normalized barrier level (r = a/x rms ), A is the probability of starting within the thresholds (A = 1 − exp(−r 2 /2)), f z is the average rate of zero crossings, and δ e is an semi-empirical measure of bandwidth (δ e = δ 1+b ). b a non-negative constant which, in this case, is equal to 0.2, and δ is a measure of bandwidth based on the spectral moments (Vanmarcke, 1972) defined as: The cumulative distribution function (CDF) of the peak values is obtained by setting t equal to D gm in equation (7); that is, the probability of the peak of the time history being less than r × x rms is equal to the probability that the time history will remain within the thresholds ±r × x rms for the entire ground-motion duration. With that, the CDF of P F is equal to: The expected value of P F can be computed with the probability density function (PDF) of P F (Equation (10)), which requires the derivation of the PDF. However, P F is continuous and defined on the positive side of the real line; thus, the expected value of P F can be computed directly from the CDF with equation Equation (11).
The mean estimate of the RVT P SA can be computed by substituting the expected value of the V 75 P F in Equation (6).

Ground-Motion Duration
In RVT, a measure of duration is needed in two steps: in the calculation of the peak factor, and in the calculation of x rms . Due to transient nature of a ground-motion, the duration measures used in these two steps are often different. D gm is the ground-motion duration, which is used in the calculation of P F ; D rms is the duration measure for the calculation of x rms . which is defined in section 3.1.4.
In seismology, the ground-motion duration is most commonly defined either as bracketed or as significant duration. Bracketed duration is the time interval between the first and last time the ground motion exceeds a threshold. Significant duration is the difference in time the normalized Arias intensity reaches two specific values. For instance, the 5 − 75% significant duration is the difference between the time the normalized Arias intensity is 5% and, the time the normalized Arias intensity is 75%. The Arias intensity is defined as integral of the squared acceleration time history: The normalized Arias intensity, also known as Husid curve, is the ratio of I a at time t over I a at the end of the ground motion: In some RVT methods, D gm is set to a measure of significant duration, but in others, D gm is treated as a free parameter with units of time. For instance, Boore (2003) used the D a0.05−0.95 significant duration as D gm , while Bora et al. (2015) and Bora et al. (2019) treated D gm as free parameter and developed a duration GMM with the goal to minimize misfit between the observed P SA and the P SA computed with RVT.
In this study, D gm is defined as an interval of significant duration. Different intervals of significant duration were tested as D gm candidates to find the one that minimized the misfit between the P SA of the used dataset (P SA N GA ) and the P SA estimated with RVT (P SA RV T ); the results of this comparison are shown in the Electronic supplement, Section S1. The D a0.05−0.85 significant duration resulted in the best fit of P SA N GA for the entire frequency range, 0.1 to 100 Hz. The Abrahamson and Silva (1996) duration GMM (AS96) was selected for estimating D a0.05−0.85 for new scenarios, as to our knowledge, AS96 is the only GMM that provides an estimate for the selected duration interval. Despite the previous results, the D a5−75 , D a5−95 , D v5−75 , and D v5−95 estimates of the Kempton and Stewart (2006) duration GMM and D a5−75 , D a5−95 , and 2D a20−80 estimates of the Afshari and Stewart (2016) duration GMM were evaluated as candidates for D gm , but the D a0.05−0.85 of AS96 resulted to a better fit of P SA N GA . The results of this comparison can be found in the Electronic supplement, Section S2.
The AS96 functional form for the mean estimate or the D 0.05−0.75 duration is: where f c is the corner frequency of the earthquake: f c = 4.9 10 6 ∆σ 10 1.5M +16.05 β is the shear-wave velocity at the source, and ∆σ is the stress drop. 1/f c is the source duration, c 1 (R rup − R c ) captures the distance dependence, and c 2 S captures the site dependence. The scaling of AS96 has a physical basis because the distance and site dependence terms are additive, instead of multiplicative, to the source duration. The rational for an additive distance dependence is that small and large magnitude earthquakes are expected to have a similar increase of duration with increasing distance due to the scattering of the seismic waves. Similarly, the duration increase due to the site effects is also expected to be independent of the earthquake size. In AS96, other interval of significant duration can be calculated with Equation (16).

Correction for non-stationarity
One of RVT's main assumptions that is violated when applied in ground motions is that the signal is stationary. Especially when predicting P SA for large T 0 , an SDOF oscillator will not abruptly stop at the end of the ground motion, instead it will have a transient decaying response, which if not considered, would lead to an overestimation of x rms . To solve this problem, Boore and Joyner (1984) (JB84) proposed to include the oscillator duration (D o ) in D rms as shown in Equation (17); D o is not included in the calculation of the P F because the response of the oscillator follows a steady decay after the end of the excitation. Liu and Pezeshk (1999) The BT15 oscillator duration model was selected for the subsequent analyses, as in preliminary evaluations, the RVT P SA estimates with BT15 provided a better fit to the recorded P SA than the alternative models. Although BT12 performed equally well in estimating the P SA of mediumto-large earthquakes, it was not selected because its is not applicable to magnitudes less than 4.

Extrapolation of EAS
To ensure that entire frequency content of the ground-motion is captured in the RVT calculations, both the ergodic and non-ergodic EAS spectra are extrapolated at low and high frequencies. At low frequencies, EAS is extrapolated to 0.01Hz with an omega-square model (Brune, 1970): where f c is the corner frequency (Equation (15)), and A f min is the amplitude of the omega-squared model at the minimum frequency of the EAS (f min ). The stress drop for the calculation of f c for the omega-squared model is estimated with the Atkinson and Boore (2011) empirical relationship. A f min is estimated based on the EAS amplitudes of 1.00f min to 1.05f min frequency bin: At high frequencies, EAS is extrapolated to 100Hz with a kappa model (Anderson and Hough, 1984): κ defines the rate of decay of the high frequencies, and A fmax is the amplitude of the kappa model at the largest EAS frequency, f max . κ can be estimated with the Ktenidou et al. (2014) κ − V S30 empirical relationship: A fmax is estimated based on the EAS amplitudes in the 0.95f max to 1.00f max frequency bin: As an example of the extrapolation procedure, the median estimate of the ergodic EAS for a M 7 event, at a R rup distance of 30km, and a V S30 value of 400m/sec is extend to high and low frequencies using the omega-squared and kappa models in Figure 5, which shows that the amplitudes of the extended frequencies are in agreement with the EAS over the usable frequency range.

RVT summary and validation
In summary, all subsequent RVT calculations are performed with: the V75 P F , the median estimate of AS96 for D a0.05−0.85 as D gm , BT15 for D rms , and the extrapolation procedure described in the previous subsection.
As a validation, Figure 6 shows the residuals between the natural-log of P SA N GA and the natural-log of P SA RV T with the recommended RV T procedure. Overall, P SA RV T is in good  Figure 7 shows the mean and the standard deviation of the residuals versus T 0 . The residuals have a positive bias at T 0 = 1 − 4sec; however, this is not propagated in the non-ergodic P SA GMM, as the GMM is developed using non-ergodic factors, which are defined in the next subsection (Section 3.2). The standard deviation or the residuals is approximately 0.2 natural-log units for the entire period range.

Non-ergodic P SA factors
The non-ergodic effects of the proposed P SA GMM are expressed in terms of a non-ergodic P SA factor (F nerg ); that is, the difference of the logs the non-ergodic P SA estimate for a scenario of interest over the ergodic P SA estimate for the same scenario (Equation (23) , which are input parameters to both the ergodic and non-ergodic EAS GMMs, but also the earthquake and site coordinates, x e and x s , which define the source, path and site non-ergodic effects in LAK21. In this formulation, F nerg captures the combined effect of all non-ergodic terms; there are no separate terms for the earthquake, path, and site non-ergodic effects.
The proposed non-ergodic P SA GMM is developed by coupling the aforementioned non-ergodic with an existing ergodic P SA GMM: where y nerg is the natural log of the non-ergodic P SA median estimate, and y erg is the natural log of the ergodic median estimate. The benefit of this approach is that it separates the non-ergodic effects from the average ground-motion scaling. F nerg does not affect the average scaling of the non-ergodic P SA GMM, as LAK21 is based on BA18, and thus, their average scaling is canceled out. Furthermore, the small bias of RV T is also canceled out in this approach, as the same RV T procedure is used to compute P SA erg and P SA nerg For the average scaling of the non-ergodic P SA GMM, y erg , we chose the Abrahamson et al. (2014) (ASK14) and Chiou and Youngs (2014) (CY14) ergodic P SA GMMs. Hereafter, the non-ergodic GMM that is based on ASK14 is called non-ergodic GMM 1 , and the non-ergodic GMM that is based on CY14 is called non-ergodic GMM 2 . The main reasons ASK14 and CY14 are selected to develop the non-ergodic GMM are: i) they were developed with the same data-set as BA18, and ii) they include complex scaling terms, such as hanging-wall effects, which can be passed to the non-ergodic GMMs. The non-ergodic P SA GMM was not developed directly with RVT and LAK21 because this approach led to an overestimation the median P SA at medium-to-large periods. Figure 8 compares the four NGAWest2 GMMs: ASK14, BSSA14, CB14, and CY14 Boore et al., 2014;Campbell and Bozorgnia, 2014;Chiou and Youngs, 2014) with the spectral acceleration response spectrum created with RVT and BA18. The NGAWest2 GMMs are in good agreement with the P SA from BA18 for the M 5 event, but the comparison worsens as the size of the earthquake increases. For periods T 0 = 2 − 4sec, for the M 8 earthquake, the P SA from BA18 is a factor of two higher than the NGAWest2 GMMs, indicating that, in this period range, BA18 has a stronger magnitude scaling than the NGAWest2 GMMs. Since LAK21 is based on BA18, a non-ergodic P SA GMM developed with RVT and LAK21 will also have a stronger magnitude scaling than the NGAWest2 GMMs. Due to the effort involved in the development of the NGAWest2 GMMs, we judge that their magnitude scaling is more likely to be correct, which is why we used the non-ergodic factors approach to develop the non-ergodic GMM; however, future studies should further investigate the cause of the different magnitude scaling.
The epistemic uncertainty of the non-ergodic P SA GMM is captured by sampling the nonergodic terms of LAK21 GMM multiple times and calculating the F nerg for each sample. As shown in the example in Section 4.1, it is important to consider the inter-frequency correlation of the non-ergodic terms, otherwise the epistemic uncertainty is underestimated.

Constant Shift and Aleatory Model
The constant shift (δc 0 ), between-event residuals (δB 0 e ), and within-event within-site residuals (δW S 0 e,s ) are estimated by fitting a mixed-effects linear model to the total residuals of the non-ergodic models: The magnitude dependence of δB 0 e and δW S 0 e,s of the two non-ergodic P SA GMMs for T 0 = 0.25sec is evaluated in Figure 9. The mean of δB 0 e and δW S 0 e,s shows no trend with M , but their empirical standard deviation decreases with M . Similarly, the R rup and V S30 dependence of the δW S 0 e,s for T 0 = 0.25sec is evaluated in Figures 10 and 11 where no significant trends are found in either the mean or the standard deviation. Figure 12 shows the estimated and smoothed δc 0 of the two non-ergodic P SA GMMs. Nonergodic GMM 2 , which is based on CY14, is only estimated up to T o = 5 sec because, at larger periods, δc 0 deviated significantly from zero.
Based on the empirical standard deviation of the non-ergodic residuals (Figure 9), both φ 0 and τ 0 are modeled as magnitude dependent (Equation (26) and (27)). Figure 13 shows the period dependence of φ 0 and τ 0 for small and large magnitudes. The magnitude dependence of φ 0 and τ 0 is more significant at small periods. The increase of the within-event aleatory variability at the small periods of small magnitudes may be caused by the radiation pattern which make the amplitude of the ground motion sensitive to the azimuthal angle. For large magnitudes, which can be thought as many small events, the radiation patterns have less impact on the ground-motion variability, because the individual radiation patterns destructively interfere with each other due to the different azimuthal angles. Similarly, the larger between-event aleatory variability at the small periods of small magnitudes is believed to be caused by differences in stress drop which shifts the ground motions at frequencies above the corner frequency of the earthquake. Due to the larger rupture dimensions of the large events, any variability in the stress drop along the rupture averages out resulting in reduced between-event variability.
The total standard deviation of the two non-ergodic GMMs are 30 to 35% smaller than the total standard deviation of the ergodic GMMs. (c) (d) Figure 13: Period dependence of aleatory model parameters. (a) period dependence of φ 0M 1 , φ 0M 1 for non-ergodic GMM 1 (b) period dependence of τ 0M 1 , φ 0M 1 for non-ergodic GMM 1 (c) period dependence of φ 0M 1 , φ 0M 1 for non-ergodic GMM 2 (d) period dependence of τ 0M 1 , φ 0M 1 for non-ergodic GMM 2 Figure 14 compares the proposed models for φ 0 and τ 0 with the standard deviations of the binned residuals for T 0 = 0.25sec. Overall, the aleatory models are in good agreement with the empirical standard deviations. The discrepancy at large magnitudes is considered acceptable, as the number of large magnitude events is small to reliably estimate the empirically standard deviation.
As a comparison with previous non-ergodic models, Figure 15 shows the total standard deviation of the two non-ergodic GMMs and the total standard deviation of the SWUS15 partially non-ergodic GMM (Abrahamson et al., 2015). The standard deviations of non-ergodic GMM 1 and GMM 2 are within the low and high branches of SWUS15 for entire period range for both small-to-moderate and large events. For small-to-moderate magnitude events and T 0 < 1sec, the total standard deviations of GMM 1 and GMM 2 are larger than the median branch of SWUS15. One possible reason for this is that σ SS of SWUS15 was estimated with magnitudes greater than 4, whereas σ 0 of GMM 1 and GMM 2 were estimated with magnitudes greater than 3 which exhibit larger variability at small periods. At large events, the total standard deviations of GMM 1 and GMM 2 are between the central and lower branch of SWUS15. The GMM 1 and GMM 2 σ 0 values are expected to be less than SWUS15 σ SS central branch because in addition to the systematic site effects, GMM 1 and GMM 2 capture the systematic source and path effects; however, the fact that the σ 0 GMM 1 and GMM 2 are larger than the lower branch of SWUS15 means that the majority of the systematic effects captured by GMM 1 and GMM 2 are related to the site effects.

Effect of EAS inter-frequency correlation in F nerg P SA
In most GMMs, the ground-motion amplitude (i.e. P SA or EAS) at every frequency is estimated independently; however, an actual ground-motion recording has peaks and troughs. That is the amplitudes of neighbouring frequencies are correlated. For instance, if amplitude of some frequency is above the average, it is likely that amplitudes of the nearby frequencies will also be above the average. This inter-frequency correlation is important in RVT, as the response of an SDOF oscillator does not only depend on the ground-motion amplitude at T 0 but also at the frequency content around T 0 . Bayless and Abrahamson (2018) showed that the P SA variability is underestimated if the inter-frequency correlation of F AS is not considered. To illustrate the effect of the inter-frequency correlation in the calculation of F nerg P SA , we applied the proposed non-ergodic GMM with and without the inter-frequency correlation in EAS. In both cases, the scenario of interest is a M 7 earthquake in Hayward Fault 8km away from a site in Berkeley, CA. The ergodic and non-ergodic EAS of the two approaches are shown in Figure 16, and the corresponding non-ergodic P SA spectra are shown in Figure 17. The non-ergodic EAS in Figure 16a are developed without inter-frequency correlation, whereas the non-ergodic EAS in figure 16b are developed using the inter-frequency correlation model in Lavrentiadis et al. (ress).
In EAS space, both approaches resulted in the same median and epistemic uncertainty range, but in P SA space, only the median is the same. The epistemic uncertainty of P SA is larger when the EAS inter-frequency correlation is considered, because if EAS is at an extreme at T 0 it will generally stay at the extreme over the neighbouring frequencies; thus, all the frequencies which influence the response of the oscillator will constructively interfere leading to a range of P SA amplitudes that is wider. In contrast, if the EAS amplitudes are uncorrelated, they will have negating effect on the response of the oscillator, resulting in a narrower range of P SA. This shows the importance of considering the EAS inter-frequency correlation in the non-ergodic P SA calculations, as otherwise, the epistemic uncertainty of the P SA is underestimated.

Magnitude dependence F nerg P SA
As an application example, Figures 18 and 19 present the EAS and P SA non-ergodic for T 0 = 0.1sec (f 0 = 10Hz) for a M 3 and M 8 earthquake in San Andreas fault. The EAS non-ergodic factors are magnitude independent; the median estimate and epistemic uncertainty of F nerg EAS is the same in both events (Figure 18). The magnitude independence allows F nergEAS to be estimated from the more frequent small magnitude earthquakes and directly applied to the large magnitude events, which are typically of more interest. This is not the case for the P SA non-ergodic factors; F nergP SA depend on the spectral shape; which is why F nergP SA are different in the M 3 and M 8 earthquakes (Figure 19), which illustrates why the non-ergodic P SA GMM is developed with non-ergodic factors that based on EAS. Most of the regional data that are used to estimate the non-ergodic effects are in form of small magnitude events, which couldn't be used if P SA non-ergodic effects were estimated directly.
In addition, Figures 18 and 19 show the spatial distribution of the epistemic uncertainty. In this example, where the location of the earthquake is fixed, the spatial distribution of the epistemic uncertainty depends on the path and site location. Both the EAS and P SA epistemic uncertainties are small near stations that have recorded past events, whereas in remote areas with no available ground-motion data to constrain the non-ergodic terms, the epistemic uncertainties are larger. The evaluation of the magnitude dependence of the EAS and P SA non-ergodic factors is further examined in Figures 20 and 21. The three scenarios in this comparison are a M 3, 5.5 and 8 event in San Andreas Fault, 105km from the site in San Francisco, CA. As mentioned previously, the non-ergodic EAS factors are the same for all three events (Figure 20), while the non-ergodic P SA factors are different, especially at small periods (Figure 21), T 0 < 0.1sec. This happens because, for f 0 > 10Hz (T 0 < 0.1sec), there is little ground-motion content in EAS to resonate the SDOF oscillator, making its response, and subsequently P SA, depended on the peak of each spectrum. Similarly, the non-ergodic P SA factors for T 0 < 0.1sec depend on the non-ergodic EAS factors at the peak of each spectrum. In this example, the M 3 event has the largest non-ergodic P SA factors at T 0 < 0.1sec, because the non-ergodic EAS factors are predominately positive over its peak (f = 2 to 6Hz). The M 8 event has the smallest non-ergodic P SA factors at T 0 < 0.1sec because its peak (f < 0.1 to 6Hz) encompasses the dip of the non-ergodic EAS factors that occur from f = 0.3 to 2Hz.

Example Hazard Calculations
A comparison of the ergodic and non-ergodic PSHA results for P SA(T 0 = 0.25sec) for a site in Berkeley, CA is presented in Figure 22. The PG&E source model was used in all hazard calculations (Pacific Gas andElectric Company (PG&E), 2015, 2017). The ergodic hazard calculations were performed with the ASK14 and CY14 GMMs, with equal weights, while the non-ergodic hazard calculations were performed with non-ergodic GMM 1 and GMM 2 , with equal weights. The epistemic uncertainty of the non-ergodic GMMs was captured by 100 realizations of F nerg P SA . This leads to a logic tree with 200 branches; each branch is a combination of a non-ergodic model (GMM 1 or GMM 2 ) and a F nerg P SA sample.
The difference between the two non-ergodic hazard calculations is that, in Figure 22b, only the regional systematic site-effects are constrained, while, in Figure 22c, recordings from past earthquakes are assumed to be available and thus, both the regional and site-specific site effects are constrained. The regional site effects are captured by the δc 1a,s term of LAK21 GMM which is a function of the site location. The site-specific site effects are captured by the δc 1b,s term of LAK21 GMM, which can be determined either from past recordings or through a site-specific site response analysis.
For the ergodic hazard calculations, the mean hazard curve is flater than the non-ergodic hazard curves due to the large aleatory variability of ASK14 and CY14, and the epistemic uncertainty is small as it only encompasses the epistemic uncertainty in the seismic source characterization and the median scaling of the ground motion averaged over all of California. Comparing the two non-ergodic calculations, the mean hazard curve is flatter and the epistemic uncertainty is larger in Figure 22b as δc 1b,s is free. This example shows the impact of non-ergodic GMM in PSHA where at moderate-to-large return periods it can lead to a factor of two to four change in the mean ground-motion level.

Conclusions
A new approach to develop non-ergodic P SA GMMs is presented in this study which considers the magnitude dependence of the non-ergodic terms. Due to the linear properties of Fourier Transform, a non-ergodic EAS GMM is used to estimate the non-ergodic effects from the small magnitude events and transfer them to the events of interest. RVT is used to compute the non-ergodic P SA effects based on the non-ergodic EAS effects, while the average scaling of the non-ergodic P SA GMM is controlled by an existing ergodic P SA GMM. Two non-ergodic P SA GMMs are developed in this study. The first one uses the ASK14 GMM as a backbone model for the average scaling and is applicable to periods T 0 = 0.01 − 10sec. The second one uses the CY14 GMM as a backbone model for the average scaling and is applicable to periods T 0 = 0.01 − 5sec. The non-ergodic P SA effects are quantified in terms of non-erodic P SA factors, that is the difference between the log of P SA estimated with RVT and the non-ergodic EAS and the log of P SA estimated with RVT and the ergodic EAS. In both cases, the LAK21 GMM is used for the non-ergodic EAS and the BA18 GMM is used for the ergodic EAS. The RVT calculations are performed with the V75 P F , the median estimate of D a5−85 from AS96 for the ground-motion duration, and the BT15 for the oscillator duration. The RVT components were chosen based on a thorough evaluation of alternative models for the peak factors, ground-motion duration and oscillator duration. The objective of the evaluation was to minimize misfit between the observed P SA and the P SA computed with RVT.
The advantages of developing the non-ergodic GMM with an ergodic backbone model and non-ergodic P SA factors, instead of developing it directly with RVT and the LAK21 are: i) the elimination of the small bias of RV T at T 0 = 1 − 4sec, ii) the separation of the non-ergodic effects from average scaling, and iii) the adoption of complex scaling terms present in ergodic P SA GMMs. Compared to the recorded P SA, the P SA estimated with RVT has a small positive bias at T 0 = 1 − 4sec. This bias is not propagated in the non-ergodic P SA factors; it is canceled out, as both the ergodic and non-ergodic RVT P SA estimates are calculated with the same approach.
Aleatory aleatory variability of the two non-ergodic P SA GMMs is approximately 30 to 35% smaller than the aleatory variability of an ergodic P SA GMM.
Future studies should reevaluate the RVT and EAS models so that when combined they result in a P SA predictions consistent with P SA GMMs. Furthermore, the proposed non-ergodic GMMs were developed with a subset of the NGAWest2 database which was compiled in 2014. As larger data sets which include more recent and more frequent small magnitude events become available, the proposed models should be assessed and potentially expanded with additional non-ergodic terms. Similarly, 3D broadband numerical simulations or inferred intensity measurements from historical earthquakes should be used to evaluate the efficacy of the proposed models.

Software and Resources
The RVT calculations were performed with the pyRVT library (Kottke, 2020) in the computer language Python (Van Rossum and Drake, 2009). The linear mixed-effects regressions were performed with the lme4 package (Bates et al., 2015) in the statistical environment R (R Core Team, 2020). The PSHA calculations were performed with HAZ45.3 (Abrahamson, 2021).

Acknowledgements
This work was partially supported by the PG&E Geosciences Department Long-Term Seismic Program. The authors thank Nicolas Kuehn for instructive comments on an early draft of this manuscript.

Declarations Funding
This work was partially funded by the PG&E Geosciences Department Long-Term Seismic Program.

S1 Comparison of different duration intervals for D gm
Figures S1 to S5 show the residuals between the records' P Sa and the P Sa estimated with RVT.
The RVT P Sa were estimated with V 75 peak factors, records' actual duration for D gm , and BT 15 for D rms . D gm is equal to: D a5−75 in Figure S1, D a5−80 in Figure S2, D a5−85 in Figure S3, D a5−90 in Figure S4, and D a5−95 in Figure S5. Figures S6 to S12 show the residuals between the records' P Sa and the P Sa estimated with RVT.