A non-ergodic ground-motion model of Fourier amplitude spectra for France

We used an ergodic ground-motion model (GMM) of California of Bayless and Abrahamson (Bull Seismol Soc Am 109(5):2088–2105, 2019) as a backbone model and incorporated the varying-coefficient model (VCM), with a modification for anisotropic path effects, to develop a new non-ergodic GMM for France based on the French RESIF data set (1996–2016). Most of the earthquakes in this database have small-to-moderate magnitudes (M2.0 – M5.2). We developed the GMM for the smoothed effective amplitude spectrum (EAS) rather than for elastic spectral acceleration because it allows the use of small-magnitude data to constrain linear effects of the path and site without the complication of the scaling being affected by differences in the response spectral shape. For the VCM, the coefficients of GMM can vary by geographical location and they are estimated using Gaussian-process regression. There is a separate set of coefficients for each source and site coordinate, including both the mean coefficients and the epistemic uncertainty in the coefficients. We further modify the anelastic attenuation term of a GMM by the cell-specific approach of Dawood and Rodriguez-Marek (Bull Seismol Soc Am. 103 (2B): 1360–1372, 2013) and Kuehn et al. (Bull Seismol Soc Am 109 (2): 575–585, 2019) to allow for azimuth-dependent attenuation for each source which reduces the standard deviation of the residuals at long distances. As an example, we compute the 5 Hz seismic hazard for two sites using the non-ergodic EAS GMM. At the 1E-4 annual frequency of exceedance hazard level, there can be a large difference between the ergodic hazard and the non-ergodic hazard if the distance between the chosen site and the available data is smaller than the correlation length of the non-ergodic terms). The combination of the non-ergodic median ground motion and the reduced aleatory variability can have large implications for seismic-hazard estimation for long return periods. For some sites, the estimated hazard will increase and for other sites, the estimated hazard will decrease compared to the traditional ergodic GMM approach. Due to the skewed distribution of the epistemic uncertainty of the hazard, more of the sites will see a decrease in the mean hazard at the 1E−4 hazard level than will see an increase as a result of using the non-ergodic GMM.


Introduction
The results of probabilistic seismic hazard analysis (PSHA) are sensitive to the standard deviation (also called sigma) for empirical ground-motion models (GMMs); even small reductions in sigma may have a significant impact on the hazard at long return periods (Bommer and Abrahamson 2006). As ground-motion data sets have grown over the past decade, there has been a trend of moving from ergodic to non-ergodic GMMs. The 2008 Next Generation Attenuation-West1 (NGA-W1) GMMs developed for crustal earthquakes (Power et al. 2008) were fully ergodic models that were applied to all regions within the same broad tectonic category. Six years later, with a much larger data set available, four of the 2014 Next Generation Attenuation-West2 (NGA-W2) GMMs for crustal earthquakes (Abrahamson et al. 2014;Boore et al. 2014;Campbell and Bozorgnia 2014;Chiou and Youngs 2014) included regional differences in four terms of the GMMs: constant term, large distance scaling, the V S30 scaling, and the basin depth scaling. Within broad regions, there are still regional differences in the median ground motion on a scale of 10 s of km. This has led to the development of non-ergodic GMMs in which the median ground motion depends on the coordinates of the source and the site.
Non-ergodic GMMs have a significantly smaller aleatory variability compared to ergodic GMMs. This reduction in the aleatory variability has been seen in many studies (Atkinson 2006;Anderson and Uchiyama 2011;Lin et al. 2011;Rodriguez-Marek et al. 2013;Kotha et al. 2016;Lanzano et al. 2017;Sung and Lee 2019) with different data sets, and have shown that for a specific site and earthquake pair, the variance of the aleatory variability is only 30-40% of the ergodic variance, indicating that most of the variability treated as randomness in ergodic GMMs is actually due to systematic source, path, and site effects. A partially non-ergodic approach (called single-station sigma) that only accounts for the systematic site terms have been used in seismic hazard studies for nuclear power plants and dams over the past ten years (e.g., Renault et al. 2010;BC Hydro 2012;Coppersmith et al. 2014;Bommer et al. 2015;Geopentech 2015;Tromans et al. 2019). A fully non-ergodic GMM for California was developed by Abrahamson et al. (2019) by combining the varying-coefficient model (VCM) of Landwehr et al. (2016) with the cell-specific approach (Kuehn et al.2019). Fully non-ergodic GMMs have not yet been used for projects, but as these models begin to be applied, there will be significant changes in the hazard for sites with nearby data (earthquakes and recording stations). For sites without nearby data, there will be no change in the mean hazard as compared to the ergodic hazard, but there will be larger epistemic uncertainty in the hazard.
Most GMMs for engineering applications are developed for 5%-damped pseudospectral acceleration (PSA); however, PSA scaling depends on the spectral shape, so the linear source, path, and site effects will not have the same scaling on the PSA values for small and large magnitudes. In contrast, the Fourier transform is a linear operator so linear effects from small magnitudes can be applied to larger magnitudes. As an example, Fig. 1a shows two sets of theoretical source spectra based on the omega-squared spectrum (Aki 1967;Brune, 1970) with the stress drop = 100 bars, kappa = 0.02 s, for M3 to M7. One set is for a hard-rock site condition with Vs = 3600 m/s (dashed lines), and one set is for the B/C boundary of Atkinson and Boore (2006) with the time-averaged shear-wave velocity over the top 30 m (V S30 ) = 760 m/s (solid lines). The site amplification function is shown in Fig. 1b. For the FAS, the linear transfer function is the same for all magnitudes; however, for the response spectrum, the site amplification depends on the magnitude due to different spectral shapes. As an example, we used the Random Vibration Theory (RVT) to convert 1 3 the FAS (Fig. 1a) to the response spectra as shown in Fig. 1c. The site amplification for each magnitude is shown as a function of frequency in Fig. 1d. This example shows that the PSA site factors depend on the spectral shape even for linear (low strain) site effects. Therefore, the non-ergodic model is developed for the Fourier amplitude spectra (FAS) values rather than for the PSA values. Developing the GMM for FAS has the added advantage that it is easier to constrain the extrapolation of the GMM using the seismological theory for ground-motion scaling (Brune 1970;Boore et al. 2014) and can provide more useful constraints for the input parameters of finite-fault simulations.
In this paper, we summarize the FAS ground-motion model developed by Sung and Abrahamson (2020) (called SA20), including the selection of the ground-motion data, the basis for the French ergodic model, and the incorporation of the non-ergodic terms to develop a fully non-ergodic GMM for France including epistemic uncertainty. In the final section, we show an example seismic hazard analysis for 5-Hz EAS and compare the nonergodic hazard with the results using a traditional ergodic GMM. Fig. 1 a Two sets of theoretical; source spectra based on the omega-squared spectrum (Aki 1967;Brune 1970) from M3 to M7. The dashed lines are for a hrad-rock site condition, and the solid lines are for a softrock site condition (B/C boundary). site amplification. (b) The site amplification function based on the B/C boundary of Akinson and Boore (2006). (c) Response spectra converted from the FAS using RVT. (d) The PSA site amplification functions for different magnitudes 1 3

Data sets
We used the processed ground-motion data from the Réseau Sismologique et géodésique Français (RESIF) data set, which includes more than 6500 recordings from 468 earthquakes recorded at 379 stations in Metropolitan France between 1996 and 2016. This dataset contains data recorded by accelerometric and broadband sensors (RESIF 1995a(RESIF , 1995b and is available as a version-numbered flat file containing the metadata, the response spectra for several damping values, and the Fourier amplitude spectra (Traversa et al. 2020). Two methods for the conversion to moment magnitude were used based on the year of the earthquake: for the earthquakes that occurred in France between 1996 to 2009, the moment magnitude was provided by the Sismicité Instrumentale de l'Hexagone (Si-Hex) project (Cara et al. 2015); for earthquakes that occurred after 2009, the moment magnitude was estimated using the conversion equation proposed by Grünthal et al. (2009aGrünthal et al. ( , 2009b) (Bremaud and Traversa 2019). Figure 2 shows a map with the event and station locations in the data set that were selected. For the ground-motion model, the site conditions are parameterized by the V S30 ; however, V S30 measurements are not available for about one-half of the recordings in the RESIF data set. For the stations without measured V S30 , we used the estimates from the global V S30 map of the United States Geological Survey (USGS) which uses the topographic slope as a proxy for the V S30 (Wald and Allen 2007; Allen and Wald 2009). As we will estimate the site-specific site terms as part of the development of the non-ergodic GMM, using proxies for the V S30 is not a key limitation of the study because the non-ergodic site terms will account for errors in estimated V S30 . The final subset consists of 6044 recordings from 463 earthquakes with the range of moment magnitude (M) between 2.0 to 5.2, rupture distance (R RUP ) between 2 and 660 km, and V S30 is between 171 and 3100 m/s (Figs. 3 and 4). The number of useable records for different moment magnitude and rupture distance ranges is shown in Fig. 3. Most of the records are in the M3.0 to M4.0 range and the 60-120 km distance bin. The total number of useable records at each frequency and also the number of recordings with and  Fig. 4. The usable frequency band was set based on the filters applied during the data processing.
The intensity measured used in the ground-motion model is the "Effective Amplitude Spectrum" (EAS) defined by Goulet et al. (2018). The EAS is an orientation-independent measure of the average horizontal-component FAS of the ground acceleration. The EAS is smoothed over a frequency band using the log 10 -scale Konno and Ohmachi (1998) smoothing window with weights defined as follows: The number of recordings with and without measured Vs30, and the final number of records after adding the global Vs30 from the USGS in which f 1 is the center frequency of the window, b is the window parameter (= 2 ∕b w ) and b w is a smoothing parameter. The Konno and Ohmachi smoothing window was selected for use in PEER projects because it led to minimal bias on the amplitudes of the smoothed EAS compared with the unsmoothed EAS. The PEER procedure uses b w = 0.0333 for the smoothing parameter because it leads to a minimal effect on the statistical moments of the EAS that are used in RVT to convert the EAS model to a response spectrum model.

Ergodic EAS ground-motion model
The ergodic EAS model for France is based on a simplified form of the Bayless and Abrahamson (2019) empirical EAS GMM for shallow crustal earthquakes in California (called BA19). The magnitude scaling in the BA19 functional form is based on the functional form used in the Chiou and Youngs (2014) GMM (called CY14) for response spectral values. The advantage of this form is that it is consistent with the magnitude scaling of FAS from small to large magnitudes (M3-M8). The functional form used for the French ergodic GMM includes the magnitude scaling, path scaling, linear site amplification, and depth to the top of the rupture (Z TOR ) scaling, and a site term as shown in Eq. 2: in which M is the moment magnitude; R rup is the shortest distance from the site to the rupture plane in km; c RB is the midpoint of the transition in distance scaling; V S30 is the timeaveraged shear-wave velocity over the top 30 m is in m/s; Z tor is the depth to the top of the rupture plane in km; B e is the between-event residual, S2S s is the between-site residuals, and WS es is the within-site residual in natural logarithm units for earthquake e and station s. The total residual, es is B e + S2S s + WS es . The total variance ( 2 ) is computed from the sum of the between-event variance ( 2 ), the within-site variance ( 2 SS ), and the betweensite variance ( 2 S2S ). The single-station total standard deviation ( SS ) sometimes used in PSHA is given by For the French ground-motion data set, the EAS is not reliable at long periods (T > 1 s), and the amplitudes are in the linear range, so the basin-depth scaling and the non-linear site effects in the BA19 functional form are not included. The maximum magnitude in the database is M5.2; that is, there are no ground-motion data from large/moderate magnitude ln EAS es =c 1 + c 2 ( − 6) + c 3 ln 1 + e c n (cM− ) + c 4 ln R rup + c 5 cosh c 6 max − c hm , 0 In the regression, we used the maximum-likelihood technique based on the randomeffects approach (Abrahamson and Youngs 1992) to estimate the coefficients for the ergodic GMM in the statistical software R (Pinheiro et al., 2020). This procedure leads to the separation of total residuals into between-event residuals, between-site residuals, and within-site residuals (Al-Atik et al., 2010). The coefficients are only estimated for frequencies greater than or equal to 1 Hz is because the number of useable recordings decays quickly for frequencies less than 1 Hz.

Non-ergodic EAS ground-motion model
Using the notation of Al-Atik et al. (2010), the median non-ergodic GMM can be written as in which erg (⋅) is the ground-motion relation of the ergodic model with a vector of predictors which includes the earthquake magnitude (M), source-to-site distance (R), and other conditions (S), t e and t s are the coordinates of the earthquake and the site, respectively, δL2L(t e ) is the adjustment for the source term based on the coordinates of earthquake source location, δS2S(t e ) is the adjustment for the site term based on the coordinates of site location, and δP2P(t e , t s ) is the path term. The hypocenter location is used for source location for these small-and moderate-magnitude events.
We relate the non-ergodic terms in Eq. 3 to the notation of Lavrentiadis et al. (2021) with modifications for the anisotropic linear R scaling: There are two constants for the site terms, one is correlated ( c stat ( ) ), and another is uncorrelated ( S ′

S2S
). There is a spatially correlated change in the V S30 scaling term based on the coordinates of site location, f 30 ( 30 ; ) = c Vs30 ln(min( 30 , 1000)∕1000)) . The ΔR Rup,i is the length of the ray path from the source to the site in the i th cell (Fig. 5). In the reference ergodic GMM, the linear R scaling coefficient, c 7 , is constrained to be negative to represent Q effects. The minus in front of c 7 R Rup term indicates the linear R scaling from the ergodic GMM is removed so that a similar constraint on the attenuation from the ergodic model ( c ca,i < 0 for each cell) can be applied to the non-ergodic model.
Here, we only model the non-ergodic effects for the linear site, path, and source because we use low amplitude ground motions. The linear site and distance terms are based on linear wave propagation through the crust that can be applied to large-magnitude events. (4) L2L = c eq ( ) The nonlinear site effects are not modeled in this approach. Any nonlinear effects need to be constrained by other GMMs that are based on observed strong ground motions or analytical modeling of site response. We also included a linear source term in the nonergodic model, but the scale factors for small-magnitude sources may or may not apply to large-magnitude sources. For example, if the average stress drop for small-magnitude earthquakes in a region is higher than average, this does not mean that stress drop for large-magnitude earthquakes in the region will also be high. In a hazard analysis, users can decide if the non-ergodic source term of the non-ergodic GMM should be applied or not. The total residuals, ε es , are modeled in two steps to estimate these non-ergodic terms: (1) the source and site terms using the VCM via a Gaussian process regression (Landwehr et al., 2016), and (2) the path term of anelastic attenuation per cell from the cell-specific approach using the Bayesian hierarchical model . These two steps are described below.

Source and site terms for the VCM
In the following sections, a brief description of the non-ergodic GMM methodology as applied for this data set is given. Readers unfamiliar with the nonergodic GMM methodology should refer to Lavrentiadis et al. (2022) for an expanded introduction to non-ergodic GMMs.
Using Eqs. 4 and 5, the total residuals from the ergodic model are modeled using the functional form in Eq. 7: in which C 0 is the intercept that accounts for the change in the implicit weighting for each recording due to the inclusion of spatial correlations of the non-ergodic terms; B e is the remaining aleatory variability of the event term; and WS es is the remaining aleatory variability of within-site term. The S ′ S2S has the same role as δS S2S for the single-station standard deviation approach, but it is computed relative to the site scaling in the VCM model.
The median ground motion using the VCM model is given by (7) es = C 0 + c eq t e + c stat t s + f Vs30 (V S30 ; t s ) + B e + S � S2S + WS es Fig. 5 Schematic showing how the length of the ray in the ith cell are calculated for the regression of the cell-specific attenuation coefficients . For the French data set, a point source (hypocenter) is used and the total residuals from the VCM are given by: The regression approach for the VCM is the Gaussian-process (GP) regression with a hierarchical Bayesian framework. The GP is a distribution over a function f(x), and its distribution is defined by a mean function and a covariance function: in which x is for the data (earthquakes and stations) and x' is the for the prediction. Because we are fitting the model to the residuals from the ergodic GMM, the mean function is set to 0 and the x = (t e ,t e ) and x' = (t e ′,t e ′) into covariance function. The covariance function model is the joint distribution of all random variables to build the distribution over a function with a continuous domain, such as the location of event or site. We set the covariance function (also known as the kernel function) of VCM to be consistent with Landwehr et al. (2016): in which 2 is the variance (amplitude) and is the correlation length (length scale); they are the hyperparameter in the kernel function, and j is the index of the coefficient. The prior distributions for hyperparameters are modeled as exponential distribution and inverse gamma distribution: We estimate the hyperparameters of the Gaussian process via the Markov chain Monte Carlo (MCMC) sampling in the program Rstan (Stan Development Team, 2020). We adopted 800 iterations per chain (4 chains in total). The mean value of the posterior distribution of the hyperparameter of Eq. 11 is listed in Table 1, from 1.0 Hz to 23.5 Hz.
The mean prediction and standard deviation of the epistemic uncertainty ( Ψ ) associated with ground-motion median predictions for different locations is given by Eq. (10) and Eq. (11) in Landwehr et al (2016): For the site terms are given by: For the source term is given by: VCM M, R, S, t e , t s , ... = erg (M, R, S, ...) + C 0 + c eq t e + c stat t s + f Vs30 (V S30 ; t s ) (9) es = ln EAS es − VCM M, R, S, t e , t s , ...
in which y is the observation dat is the mean prediction, Ψ is the standard deviation of epistemic uncertainty for site locations or source locations. x * denotes the new predictor variable for a new location t * . x i and x l are the ith and lth predictor variable for the existing location t i and t l for sources or sites, respectively. For example, if we apply Eq. 14 to estimate the site terms, the location t i and t l should consider using the site locations, in contrary, if we calculate the source terms, these locations should be the source locations. Moreover, n is the number of data points, so if we calculate the source terms, the n should be the number of earthquakes. κ m is from the Eq. 11, I is the identity matrix, d is the number of non-ergodic coefficients in the model. Here, the y e and y s are the between-event and with-event residuals, 0 and σ 0 are the standard deviation of aleatory variability of the y e and y s . In our study, the spatially varying coefficients ( c eq , c stat and c Vs30 ) only have epistemic uncertainty; it will not like the ground motions contain the epistemic uncertainty and aleatory variability both. For the median model, we only need the non-ergodic terms, so the 0 and σ 0 of the ground motion is not included in calculating the epistemic uncertainty in Eqs. 14-17. The non-ergodic coefficients at the locations (e.g., t i and t l ,) are not known, so we assumed the non-ergodic coefficients follow a multivariate normal posterior distribution: in which c j,n is the non-ergodic coefficient at the nth existing location, c j,n is the posterior mean, and Ψ c j,n is the posterior variances. The epistemic uncertainty of c j,n can be accounted in predicting c * j by using the marginal distribution of c * j . Based on these assumptions, the mean prediction and the epistemic uncertainty associated with nonergodic coefficient predictions for new locations is given by (Lavrentiadis et al., 2021) (16) c eq = k T K + 2 0 I −1 y e , We adopted the Eqs. 22 and 23 to calculate the mean prediction and epistemic uncertainties of the c eq , c stat , and c Vs30 for each source location and site location are estimated. For the non-ergodic PSHA calculations, these epistemic uncertainties of the adjustments for each source cell and site location are considered in the logic tree with a map of the non-ergodic terms on each branch of the logic tree for the GMM.
To show the adjustment map of source constant ( c eq ), site constant ( c stat ), and ln(Vs 30 ) coefficient ( c Vs30 ), we divide France into cells of size 0.2 × 0.2 degrees. The correlation lengths shown in Table 1 are larger than 0.2, so this grid size is adequate.
The mean values of the three non-ergodic coefficients are computed for each cell using either the source location or site location. Maps of the spatially varying mean coefficients are shown in Fig. 6 for 1 Hz and 5 Hz. These mean adjustment terms are zero in regions without data (or with sparse data) and potentially show the positive values or negative values for cells close to observed events or seismic stations.

Path-specific attenuation term
The same 0.2 × 0.2 degrees grid size is used for the cells for the non-ergodic path effects. Figure 7 shows the number of rays that pass through each cell at 5 Hz. For eastern and southern France, there are data to provide good constraints on the path effect per cell. We apply the cell-specific anelastic attenuation terms to replace an ergodic term in the GMM and use the same regression model (Bayesian hierarchical model) and settings of the prior distributions following Kuehn et al. (2019) in program Rstan. The setting are the betweenevent term, between-site term and within-site term of a non-ergodic GMM. The c 7 R rup term is the linear R scaling of the ergodic GMM (Eq. 2). To avoid nonphysical attenuation coefficients, the c ca,i terms are constrained to be negative. Figure 7 also shows the mean of the posterior distribution of the c ca,i per cell at 5.0 Hz. There are larger mean anelastic attenuation coefficients (lower Q) in the Alps region and smaller mean values (higher Q) in western France. If a cell does not include any ray paths, the mean cell-specific anelastic attenuation coefficient will become equal to the ergodic anelastic attenuation value (c 7 in Eq. 2). Sung and Abrahamson (2020) show a comparison of the Q models from Campillo and Planter (1991) and Mayor et al. (2018). These Q studies show lower Q values in the south-eastern regions and higher Q values in the western areas consistent with the non-ergodic terms. The standard deviation of the epistemic uncertainty of the c ca,i is captured from the posterior distribution. As with the epistemic uncertainties of VCM, these epistemic uncertainties for each source cell are included in the logic tree in the non-ergodic PSHA calculations.

Residuals and standard deviation
The residuals of the ergodic and fully non-ergodic (combine VCM and path effects) GMMs for 5.0 Hz are shown in Fig. 8. The results show a significant reduction of the variability of the between-event and between-site residuals for the non-ergodic model as compared to the ergodic model. The main reason for an inflated estimate of aleatory variability of the ergodic GMM is that systematic source, site, and path effects are assumed to be random and (24) apply to all sites. In contrast, the non-ergodic model incorporates the repeatable systematic source, path, and site effects into the GMM. At short distances, the variability of the within-site residuals is similar for the two models, whereas there is a significant decrease in the variability at long distances for the non-ergodic GMM due to the path terms. Figure 9 shows the aleatory standard-deviation terms for the ergodic, partially non-ergodic, and fully non-ergodic GMMs for each frequency. There is a significant Fig. 8 Residual of the ergodic GMM (triangle) and non-ergodic GMM (circle) for 5.0 Hz reduction of the aleatory variability for the non-ergodic model. For example, at 5 Hz, the ergodic aleatory standard deviation is 0.94, the partially non-ergodic aleatory standard deviation (the single-station standard deviation) is 0.76, and the non-ergodic aleatory standard deviation is 0.59. This corresponds to a 60% reduction in variance from the ergodic to the non-ergodic GMM and is consistent with the results from previous studies for active regions discussed in the introduction.

Median EAS prediction and uncertainty for non-ergodic model
Because the adjustment terms of a non-ergodic GMM vary by the geographical location, the median predictions at a given site also change depending on event locations. Figure 10 shows an example of the spatial variation of median EAS predictions by the French nonergodic GMMs for a set of predictor variables M = 5.0, V S30 = 2100 m/s, Z TOR = 10 km at 5.0 Hz for Site1 in southeastern France (43.6748°N, 5.7664°E). Comparing these two Fig. 9 The aleatory standard deviation (ln unit) of ergodic, partial, and non-ergodic GMMs models, the ergodic GMM shows the same attenuation in all directions, but for the nonergodic GMM, the attenuation can depend on the direction based on the non-ergodic terms. Figure 11 compares the distance scaling of the EAS in different directions for the ergodic and non-ergodic models at 5.0 Hz for two sites: Site1 and Site2 (47.2294°N, 0.1673°W). The largest differences are seen for distances of 100-300 km due to the path effects. The standard deviation of the epistemic uncertainty of the median ground motion for the non-ergodic model for two sites is shown in Fig. 12. This epistemic uncertainty is from the combination of the standard deviation of the posterior distribution for all the nonergodic terms (source, site, and path). The results show the lower uncertainty is constrained to locations where data are available, whereas the uncertainty is larger for the region with sparse or no data. All epistemic uncertainties for Site 1 are lower than the resulting of Site 2, which reflects the smaller uncertainty of the site term for Site 1 than for Site 2.
The epistemic uncertainty versus the number of recordings of Site 1 is shown in Fig. 13 for all the source grid and the source grid only within the source-to-site distance of 0-km to150-km. The epistemic uncertainty in the non-ergodic terms can reflect the number of recordings near the site of interest. That is, for the area without data, Fig. 11 For 5 Hz-results, a northern/western source grids for Site 1 and b southern/eastern source grids for Site 2. c The distance scaling of the ergodic/non-ergodic EAS predictions (g.s) for Site 1. d The distance scaling of the ergodic/non-ergodic EAS predictions (g.s) for Site 2 1 3 the non-ergodic GMM can still be applied, but the epistemic uncertainty is large. For example, in Site 1, the epistemic uncertainty is about 0.425 when the source grid did not include the available observation. In contrast, the epistemic uncertainty can decrease to 0.3-0.32 when the source grid with 3-20 recordings. We use the same parameters (V S30 = 2100 m/s, Z TOR = 10 km) to calculate the median EAS spectra at 1.0 Hz to 23.5 Hz for two target sources which have the same source-site distance of 247 km. Figure 14a shows the first target source (Source 1, 43.7°N, 2.7°E) is located west of the site, the second source (Source 2, 45.9°N, 5.7°E) is located north of the site, and the third source is located northeast of the site (Source 3, 44.1°N, 6.1°E). Figure 14b compares the median EAS spectra from the non-ergodic and ergodic models from M3 to M7 for source 1 and source2, which have similar source-site distances. For Source 1, the ergodic and non-ergodic models have similar Fourier spectra, but for Source 2, there is a significant difference at frequencies greater than 2 Hz between non-ergodic and ergodic models. Source 2 has observed data nearby, so it shows a large adjustment for the non-ergodic model. Figure 14c only shows the source 3 results that is nearby Site1 (the distance of 40 km). In this target source, there is an apparent adjustment between the ergodic model and non-ergodic version for all the frequency dues to nearby observations. Figure 14d shows the epistemic uncertainties of the non-ergodic model at all frequencies for each source. Source 2 and Source 3 have smaller epistemic uncertainties than Source 1 due to the available earthquakes near these two sources.

Examples for non-ergodic hazard calculations
For PSHA implementation of the non-ergodic approach, the site/source specific adjustment along with the standard deviation of epistemic uncertainty in the adjustment needs to be included. The current approach is to precompute the net adjustment term (sum of the non-ergodic terms) for each source location for a single site location. We randomly sample 100 different adjustments maps from the distribution of the epistemic uncertainty of the non-ergodic terms. Figure 15 shows the two samples of 100 maps for sampled adjustments, the range of latitude and longitude is from 40°N to 52°N and 6°W Median EAS (g.s) spectra of ergodic and non-ergodic models from M3 to M7, for Source 1 and Source2, and c for Source 3. d Epistemic uncertainty (ln unit) per frequency for Source 1, Source 2 and Source 3 to 12°E, respectively. In PSHA calculation, we add these adjustments to the median ground motion from the ergodic model. In this approach, there is a logic tree with 100 branches of the non-ergodic terms for the GMM that are equally weighted (Fig. 16).
For each realization, a single non-ergodic term is estimated along with a set of spatially correlated non-ergodic source and path terms. The net adjustment is the sum of the source, path, and site terms. A total of 100 realizations are generated to sample the variance and sptial correlation of the non-ergodic terms. By precomputing the net adjustment terms, fewer changes are required for the PSHA program. The main change is that the adjustment terms depend on the latitude and longitude of the source.
The following is the formula for the ergodic probability of exceeding any ground motion level using knowledge of the median ground motion ( erg ) and ergodic aleatory standard deviation ( erg ) in the units: where the Y is the ground-motion parameter of interest, and P(Y > Z|M, R) is the conditional probability that Y is larger than Z for a given magnitude (M), distance (R), and other relevant parameters (e.g., T (period)). Φ() is the cumulative distribution function of a standard normal distribution. But for non-ergodic GMMs, the latitude and longitude of each source need to be passed to the GMM subroutine and the appropriate adjustment factor interpolated from the pre-calculated values.
in which Lat e and Lon e are the latitude and longitude for the source, and Lat s and Lon s are the latitude and longitude for the site. The HAZ45 code runs the hazard for a single site, so the Δ ne is a table of adjustments for different source locations and magnitudes for the selected site, and the ne is the non-ergodic aleatory standard deviation. Multiple realizations are included to capture the epistemic uncertainty in the non-ergodic terms. Site 1 and Site 2 per source grid. For sites without local data (e.g., Site 2), the mean values in the maps are close to zero. For sites with location data to constrain the non-ergodic terms (e.g., Site 1), the mean can be close to zero or different from zero depending on the data. The epistemic uncertainties range from 0.3 for regions with local data to 0.7 for regions without local data. It should be similar to the previous epistemic uncertainty map (Fig. 12) because we adopted these epistemic uncertainties to generate the samples (adjustments) for each cell. Figure 18 compares the hazard using the ergodic EAS GMM with the hazard from the 100 samples of the non-ergodic adjustment terms for two sites, the seismic source model is the EDF zoning (Drouet et al., 2020): Site1 (close to the local data) and Site2. For Site 1, there is a large reduction in the mean non-ergodic hazards for 5 Hz. At an annual frequency of exceedance of 1E-4, there is a factor of 1.6 reduction in the ground motion, and the 95th fractile from the non-ergodic hazard is close to the mean ergodic hazard. The slope of the median hazard is steeper than the slope of the mean hazard. For Site 2, the 5-Hz mean hazard are similar to the mean ergodic hazard results. The adjustments for Site2 are small due to the sparse available data in this region, and there is increased epistemic uncertainty leading to a wider range of the fractiles than Site1.
For the non-ergodic GMM, if only the aleatory variability is reduced, then the nonergodic hazard curves would be much steeper than the ergodic hazard curves; however, the adjustments vary spatially, which acts like additional variability in the total hazard (summed over all source locations). So, the non-ergodic hazard curves are not as steep as would be the case if the adjustment was the same for all source locations, as is the case with the partially non-ergodic single-station sigma approach. The mean hazard also samples the epistemic uncertainty from the 100 realizations which will also flatten the slope of the mean hazard curve for large epistemic uncertainty. The effect of the epistemic uncertainty on the slope of the hazard curve can be seen by comparing the slope from the mean hazard and the slope from the median hazard.

Fig. 18
Mean hazard as well as 5%, 50%, and 95% fractiles of the resulting hazard curve distribution for Site1 (Left) and Site2 (Right) for 5.0 Hz

Conclusions
Traditionally, GMMs used in PSHA have been based on the ergodic assumption (Anderson and Brune 1999). The move to non-ergodic GMMs will lead to GMMs that have a better physical basis and provide more accurate estimates of the hazard in regions with local data. In regions with little local data, non-ergodic GMMs can still be used. In this case, the mean hazard will approximate the ergodic mean hazard, but the uncertainty range will be broader and will provide a better estimate of the uncertainty and the value of collecting local ground-motion data.
The non-ergodic GMM was developed for EAS rather than PSA to avoid the scaling issue being affected by differences in the response spectral shape which allows the use of non-ergodic terms estimated from small-magnitude earthquakes to be applied to large-magnitude earthquakes. The non-ergodic and ergodic GMMs for EAS developed in this study can be converted to GMMs for response spectral values using RVT. An example of using traditional RVT to convert an EAS GMM to a PSA GMM is given by Lavrentiadis and Abrahamson (2022). An empirically calibrated RVT approach is given by Phung and Abrahamson (2022). In these two approaches, RVT is applied to both the ergodic EAS GMM and the non-ergodic EAS GMM. The non-ergodic/ergodic PSA ratio is then used to adjust the ergodic PSA GMM to be a non-ergodic PSA GMM. This approach has the advantage that systematic biases in the RVT method cancel by using the ratio, leading to more stable results. For the French model, the extrapolation to low frequencies should be revised befor applying RVT. Because the available data groundmotion from France were not reliable at low frequencies, the current model assumed that the low-frequency scaling from the BA19 model is applicable to France at frequencies less than 1 Hz, which may not be appropriate for France.
Using a non-ergodic GMM, there is a large reduction in the aleatory standard deviation, but there is also a shift in the median from each source location that can be either positive (large median) or negative (smaller median). At the 1E-4 annual frequency of exceedance, there tends to be a reduction in the hazard using the non-ergodic GMM, but at some sites, there will be an increase in the hazard if there is an increase in the medians for the controlling sources that offsets the reduction in the aleatory standard deviation at the 1E-4 annual frequency of exceedance.