Lessons learned from applying varying coefficient model to controlled simulation datasets

The development of site- and path-specific (i.e., non-ergodic) ground-motion models (GMMs) can drastically improve the accuracy of probabilistic seismic hazard analyses (PSHAs). The varying coefficient model (VCM) is a novel technique for developing non-ergodic GMMs, which puts epistemic uncertainty into spatially varying coefficients. The coefficients at nearby locations are correlated by a prior distribution imposed on a Gaussian Process. The correlation structure is determined by the data, and later used to predict coefficients and their epistemic uncertainties at new locations. It is important to carefully verify the technique before its results can be accepted by the engineering community. This study used a series of simulation-based controlled ground-motion datasets from CyberShake to test a modified VCM technique, which partitions the epistemic uncertainty into spatially varying source, site, and path terms. Because the simulation parameters (inputs) are known, verification of what is recovered by the VCM from CyberShake simulations is straightforward. We found that the site effects in CyberShake datasets can be reliably estimated by the VCM. However, the densely-located self-similar events in CyberShake datasets along pre-defined faults violate the isotropic assumption underlying the VCM, thus preventing the VCM from capturing the genuine source effects. For path effects, cell-specific attenuation approaches fail to recover the anelastic attenuation pattern of the 3D velocity model, which is most likely due to other unmodeled effects and inappropriate assumption of wave-propagation path. Instead, a midpoint approach that only considers the aggregated path effects can better recover the strong anelastic attenuation within basins by fixing the correlation length of path effects. Lessons learned in this study not only provide guidance for future applications of VCM to both simulation and empirical datasets, but will also guide further development of the technique, with emphasis on the recovery of path effects.


Introduction
and perform preliminary applications on various datasets; this entire special issue captures the work completed by the NEWG in the last few years, and introduced by Lavrentiadis (2021). In coordination with the NEWG, our role is to perform verification and sensitivity studies of the VCM technique using simulated ground-motion datasets. There are several major limitations in current ground-motion simulations. For example, the epistemic uncertainty in wave-propagation through 3D velocity models is large and difficult to constrain. The earthquake-source effects in simulations tend to be simpler than for real earthquakes, but it is prescribed and therefore easier to define numerically. Moreover, simulated groundmotions are not constrained to be applicable to the Gaussian process, thus applying the VCM may lead to biased results. On the other hand, ground-motion simulation datasets present two major advantages over empirical datasets. First, simulations can be obtained for any source regions, sites, and paths, all of which can be sampled many times. Hence, compared to the sparse empirical datasets, the source, site, and path effects can be more accurately estimated. Second, since the source parameters of earthquakes, site conditions, and the underlying velocity model used in simulations are known, the overall accuracy of the VCM outputs can be directly verified against the inputs.
This study used the Southern California Earthquake Center CyberShake ground-motion simulation platform (hereafter CyberShake), which is, by definition, a non-ergodic PSHA model that incorporates an ERF with a full 3D wave-propagation to perform physicsbased PSHA by computing ground-motions at chosen sites Jordan et al. 2018). At each site, seismograms are calculated using seismic reciprocity for all ruptures within 200 km from the Uniform California Earthquake Rupture Forecast, Version 2 (UCERF2) model (Field et al. 2009). The ruptures themselves are generated using a kinematic source model Pitarka 2010, 2015) and the wave-propagation properties are based on large regional 3D community velocity models (CVMs). CyberShake includes up to 415 variations (i.e., rupture scenario) for each rupture (with the same magnitude and rupture plane) defined in UCERF2. Each rupture variation is essentially an unique event in the CyberShake datasets. Once the ruptures and their variations are defined, CyberShake uses anelastic wave-propagation simulation codes to calculate strain Green tensors (SGTs) around the site (Cui et al. 2013). Seismic reciprocity is used to post-process the SGTs and obtain synthetic seismograms. The wave-propagation and site effects are therefore assumed to be linear. RotD50 pseudo-spectral acceleration (PSA) at periods of 2, 3, 5, and 10 s are calculated from each simulated seismogram. We adopted a modified VCM technique largely based on Landwehr et al. (2016). Instead of directly solving for a non-ergodic GMM with spatially-varying coefficients, we first obtained an ergodic GMM through linear regression. Following the regression, we applied the VCM to the total residuals and computed the spatially varying source, site, and path effects. The advantage of the modification is two-fold (1) it significantly reduces the computational cost by decreasing the number of spatially varying components from a dozen to three, which is a primary concern for very large simulation datasets; and (2) the outputs have a similar decomposition to those from Al Atik et al. (2010), thus the results can be easily understood by the earthquake engineering community and compared to previous studies (e.g., Villani and Abrahamson 2015;Meng et al. 2022).
The modified VCM was applied to a CyberShake study in Southern California (hereafter CS15.4). CS15.4 includes 336 sites, 360,472 events and 97,214,974 seismograms (Fig. 1a), from which we selected subsets of simulation data. It would be impractical to apply the VCM to the entire CS15.4 dataset, as it requires an extraordinary amount of memory and computation time. Moreover, the main purpose of this work is to verify the VCM and test its sensitivity with a controlled dataset, and not to develop a non-ergodic GMM based on CyberShake data. Therefore, we first randomly generated four reduced CS15.4 datasets of different sizes. Then, the VCM was applied to the four datasets, and the outputs were compared among themselves to check the sensitivity on dataset sizes. Next, we identified results that were consistent with the input components, which can be used to better constrain future PSHA applications. When discrepancies emerged between inputs and outputs, we explored their causes by conducting additional tests on controlled datasets and proposed potential solutions for future VCM applications.

Methodology
We defined a simple function form to develop the ergodic GMM, which includes magnitude (M) scaling, geometric spreading, magnitude-distance interaction and linear site response (e.g., Chiou et al. 2008): where h is set to 4.5 km and c i are coefficients determined through regression; R rup is the closest distance between the site and the rupture plane; f(V S30 ) is the linear site response term based on the slowness-averaged shear-wave velocity in the top 30 m of the velocity model (V S30 ): The separation point of f(V S30 ) is set as 1500 m/s as it is the upper limit of the V S30 values from the NGA-West2 dataset (Ancheta et al. 2014). While we use the V S30 term here, Cyber-Shake does not model the near-surface materials to allow the simulations to capture fine-scale variations in site conditions. Hence, we use the value of the top 50 m (half a mesh size) as the V S30 value for computation purposes.
The total residuals are decomposed into three spatially varying components (source, site, and path effects) and the remaining aleatory residuals: (1) where L2L(�� ⃗ x e ) are the source effects that spatially vary with event locations �� ⃗ x e , S2S(�� ⃗ x s ) are the site effects that spatially vary with site locations �� ⃗ x s , P2P(�� ⃗ x e , �� ⃗ x s ) are the spatially varying path effects; B 0 e and WS 0 e,s are the remaining residuals after all spatial correlations are considered.
For this paper, we experimented with three modeling approaches for P2P(�� ⃗ x e , �� ⃗ x s ) , all of them being simplifications relative to the complex source-site travel paths. The first approach is the two-dimensional (2D) cell-specific attenuation (hereafter 2D cell approach), similar to the one introduced by Dawood and Rodriguez-Marek (2013). The entire study region is divided into 2D cells of 0.25 ० × 0.25 ० . The path effects are the sum of anelastic attenuation from all cells along a horizontal straight line from the closest point to the site: where δ i is the attenuation per kilometer of travel in the i th cell, which is intended to capture the anelastic attenuation in traditional GMM development. δ i is computed as an independent random effect during regression (i.e., no spatial correlation). ΔR i is the approximate travel distance within the i th cell. The second approach is the three-dimensional (3D) cell-specific attenuation (hereafter 3D cell approach). The entire subsurface is divided vertically into 1 km layers. In the top 10 km, each layer is further divided into cells of 0.25 ० × 0.25 ० . Below 10 km, each layer is no longer divided into small cells, because there is very little heterogeneity at larger depths in the 3D velocity structure used for the simulations. This implementation provides the additional advantage of significantly reducing the computational demands. The wave-propagation is assumed as a straight line from the hypocenter to the site. The third approach only considers the aggregate path effects along the entire path. Similar to the path term described in Al Atik et al. (2010), each path between a source region and a site is considered unique. We used the midpoint t p between the site and the closest point of the causative fault to represent the approximate location of a path: where 1 (����� ⃗ x mp ) is anelastic attenuation that spatially varies with midpoint locations ����� ⃗ x mp . Hereafter this approach will be referred to as the midpoint approach.
The next step is to determine the amount of correlation in L2L(�� ⃗ x e ) , S2S(�� ⃗ x s ) and . This is achieved by imposing a prior distribution for the Gaussian Process (Rasmussen et al. 2006) prior on the three terms. We used the Matern kernel to constrain the spatial correlation between any two data points: where Γ( ) is the gamma function, K is the modified Bessel function of the second kind; r is the distance in degree between two data points; is the smoothness parameter, which is fixed at 1 in this study; l and are hyperparameters that represent the correlation length and standard deviation, respectively, which are determined by the data. The Markov-chain Monte Carlo (MCMC) approach is often used to solve the hyperparameters of a Gaussian Process, but it tends to be very computationally expensive because it computes the joint posterior distribution of the hyperparameters. The integrated nested laplace In R-INLA, the spatial correlation of data is estimated by stochastic partial differentiation equation (SPDE) (Lindgren et al. 2011), which uses a finite element method to represent a Gaussian Process. First, a mesh consisting of triangles is created over the study region. The continuous spatial processes are then discretized onto the vertices of the triangles. If a data point is located within a triangle, three vertices are directly related to the data point weighted by their respective distance to it. If a data point is located on the edge of a triangle, only two vertices at two ends are related to it and weighted by distances. If a data point is located on one vertex of a triangle, the other two vertices are not related to it. The total weight of each data point is always 1. Next, the SPDE algorithm calculates the correlation structure among all the vertices of the mesh, which is then used to make predictions at any location. Smaller triangles lead to higher precision but require larger memory and longer computation times. After many trial runs with the reduced CyberShake datasets, the most satisfying triangle size for the balance of the precision and computational cost was selected to be 0.1 ० .

Datasets
The magnitude-R rup and V S30 distributions of CS15.4 are shown in Fig. 1b, c, respectively. While spatially distributed small magnitude earthquakes dominate empirical datasets, they do not exist in CyberShake. The majority of events in CS15.4 have magnitudes above 6 and all occur along pre-defined faults. To test the sensitivity of VCM to dataset size, we generated four reduced datasets by randomly selecting 100, 200, 600, and 1000 events in CS15.4 (Table 1 and Fig. 2). As a result, there are almost no overlapping events between any two reduced datasets. While these are reduced datasets relative to the whole CS15.4 simulation set, the corresponding number of data is still quite large (Table 1); in comparison, the NGA-West2 database contains 21,336 records spanning 599 events. Since a single rupture in CyberShake includes up to hundreds of variations, a single rupture variation per rupture was selected so that the reduced datasets have the largest variability in source location and parameters. There are 336 sites in CS15.4, which are all included in the reduced datasets. The sites are concentrated within a box centered around Los Angeles, which covers the Los Angeles and Ventura sedimentary basins and their immediate vicinity; two additional test sites are located away from the box. The 3D velocity model used for the CS15.4 simulations is CVM-S4.26-M01, obtained from full 3D waveform tomography (F3DT) (Lee et al. 2011;Small et al. 2017), supplemented by borehole information where available, and using a tapering scheme (Ely et al. 2010) to create a geotechnical layer that reflects the near-surface velocities from Wills (2006). The shear-wave velocity (V S ) has a floor value of 500 m/s. Outside the range of CVM-S4.26-M01, a generic 1D velocity model is used. We used RotD50 PSA at 5 s as the ground motion measure, the same period used for the development and validation of the 3D velocity model (Lee et al. 2011;Lee and Chen 2016).
The effect of anelastic attenuation is usually specified with quality factor, Q. However, because CVM-S4.26-M01 only provides the seismic velocities and density of the material, we used an empirical relationship to invert Q from the seismic velocity model (Taborda and Bielak 2013): The Q S maps are then used as the reference point for the input attenuation (Fig. 3). Near the surface, there are two regions with much stronger attenuation than the rest of Southern California. The first region includes the Ventura and Los Angeles basins, as well as their immediate vicinity; the second one is the Imperial Valley to the south. The areas of the two regions quickly decrease with depth, such that below 5 km the attenuation within the two regions becomes very similar to most of the study region, leaving the Great Valley (to the north) as the only region with stronger attenuation.

Results
There is no apparent trend in the total residuals relative to any of the predictor variables for the four datasets ( Figure  without events (Fig. 4b, d). Moreover, there is no longer a clear distinction in L2L (�� ⃗ x e ) between locations with and without data for the three larger datasets (Fig. 5b, d). In other words, the VCM fails to capture the genuine source effects when the number of events in CyberShake datasets increases, and the resultant L2L(�� ⃗ x e ) and L2L (�� ⃗ x e ) are no longer meaningful.
The results of S2S(�� ⃗ x s ) and its predictive uncertainty S2S (�� ⃗ x s ) are almost identical across four datasets (Figs. 6, 7). Within the Ventura and Los Angeles basins, the sites have the strongest site responses. In comparison, the hard-rock sites at the San Gabriel Mountains have much weaker site responses. The site responses return to zero outside of the site coverage, where they show increased epistemic uncertainty. The few test sites located outside of the box show the same trends. The consistent pattern suggests that the genuine site effects of CyberShake datasets are accurately recovered by the VCM technique.
For the 2D cell approach, the variations of δ i are constrained at locations with data (Fig. 8), and the uncertainty of δ i increases at locations without data (Fig. 9). However, the approach failed to capture the strong anelastic attenuations expected within the Ventura and Los Angeles basins observed in the Q S maps (Figs. 3, 8). In all four datasets, cells within the basins and their immediate vicinity (hereafter inner cells) have δ i close to zero (Fig. 8), which represents the median attenuation of all the cells sampled by seismic waves. Outside the basins (hereafter outer cells), the δ i patterns are dependent on the data distribution. The δ i results obtained by the 3D cell approach show a slight improvement over the 2D cell approach. At shallow depths (0-3 km), the outer cells have consistently weak attenuation, especially for the three larger datasets, which agree with the Q S maps. However, δ i of the inner cells at shallow depths are still close to zero for four datasets (Fig. 10). The number of inner cells with close-to-zero δ i quickly decreases with depth. Below 5 km, there is no longer a distinct area with close-to-zero δ i for all the datasets. Cells with weak and strong attenuation are located randomly, with very little agreement with the Q S maps at the corresponding depth.
In comparison with the 2D and 3D cell approaches, the midpoint approach enables the recovery of attenuation within the basins for the four datasets (Fig. 11), and the predictive uncertainties P2P (����� ⃗ x mp ) increase at locations without data (Fig. 12). However, similar to L2L(�� ⃗ x e ) , the variations of 1 (����� ⃗ x mp ) extend to large areas without data due to the very large correlation lengths (Fig. 11). Moreover, the correlation lengths are evidently different among the four datasets, which result in significantly different patterns outside the basins. Therefore, although the midpoint approach leads to results that are more consistent with the input velocity model, issues on the stability of the approach for different sizes remain.

Discussion
This section discusses the results for each of the spatially variable results in more detail. One point of note is that in addition to these interpretations, it is possible that the approach we used to perform a linear regression first and then apply the VCM approach to the residuals may introduce a bias. However, because the CS15.4 is large and complete (all events have estimates at all sites and vice-versa), we do not believe this to be a source of concern.

Source effects ıL2L(�� ⃗ x e )
The very large correlation lengths of L2L(�� ⃗ x e ) are likely caused by two unique properties of the CyberShake datasets. First, the majority of events are self-similar, that is, they follow the same magnitude-rupture area relationship (i.e., same stress drop) (Somerville et al. 1999). Second, the majority of events in CyberShake datasets are associated with principal faults in California. When the number of events becomes large, they occur densely along those faults and form large linear features (Fig. 2). The self-similar events along large linear features thus lead to large correlation lengths of L2L(�� ⃗ x e ) . In R-INLA, the Gaussian Process is assumed isotropic, that is, the process only depends on the distances between To test this hypothesis, we randomly selected another 600 events in CS15.4 (hereafter, CS15.4-600b) under two restrictions: (1) the minimum hypocentral distance between any two events is 10 km; and (2) the stress drops are evenly distributed between 0 and 3 MPa (Fig. 13a). After applying VCM to CS15.4-600b, we find that the correlation length of L2L(�� ⃗ x e ) becomes much smaller (Fig. 13b). The variations of L2L(�� ⃗ x e ) are constrained at event locations and correlate well with stress drops (Fig. 13b). As intended, L2L (�� ⃗ x e ) also increases at locations without data (Fig. 13c). This test confirms that large correlation lengths of L2L(�� ⃗ x e ) for CyberShake datasets are indeed caused by the dense linear distribution and uniform stress drop of events. Thus, for future application of VCM on CyberShake datasets, it is recommended to fix the correlation length at a very small value informed by semivariogram model (e.g., Walling and Abrahamson 2012), so that the source effects are constrained at event locations. This is essentially the same approach that computes event terms as point estimate random effects; see Al Atik et al, (2010).
In typical empirical datasets, earthquakes tend to be much more sparse and have a wide range of stress drops. It is also worth noting that the magnitudes in empirical datasets tend to be much smaller than for the CyberShake datasets. Therefore, it is unlikely to run into the same problem of large correlation lengths. For example, the source effects obtained by VCM for California and France datasets do not exhibit such behaviors (Lavrentiadis et al. 2021;Sung et al. 2021). However, if an empirical dataset exclusively consists of one earthquake sequence or swarm, in which all the events have similar stress drop and are located densely along one fault, the VCM will likely yield large correlation lengths of L2L(�� ⃗ x e ) and unreasonable predictions at locations without data.

Site effects ıS2S(�� ⃗ x s )
The large site responses observed within the Ventura and Los Angeles basins agree with the strong basin amplification observed from empirical data (e.g., Filippitzis et al. 2021;Wald and Graves 1998). VCM does an excellent job capturing the basin effects from the Cyber-Shake datasets. Provided that the CyberShake site effects are validated, they could be used directly to complement empirical GMMs and improve seismic hazard analysis within the site coverage. Based on the results of L2L(�� ⃗ x e ) , it is reasonable to assume that if most of the sites are in a dense linear array, the isotropic assumption of the Gaussian Process will likely yield a large correlation length of S2S(�� ⃗ x s ) as well. Hence, for future non-ergodic GMM development, one should take extra caution when dense linear arrays are involved.

Path effects ıP2P(�� ⃗ x e , �� ⃗ x s )
Both the 2D and 3D cell approaches failed to capture the strongest attenuation at shallow depths within the Ventura and Los Angeles basins for any reduced dataset. Although the cell-specific attenuation term is designed to recover anelastic attenuation, in reality it includes all the systematic effects that are not modeled in the GMM. First of all, the complicated wavefields due to the 3D velocity model are smoothed (e.g., strong reflections and refractions at interfaces where seismic wave velocities change significantly, and generation of large surface waves). Second, due to the dense site distributions in CS15.4, the azimuths from the source location to all the sites are very similar for many events. Therefore, some source effects (e.g., radiation patterns) are likely to be incorrectly mapped into the path effects. Third, the epistemic uncertainty of the 3D velocity model itself is an important component of the path effects but not modeled yet.
To test this hypothesis, we designed controlled datasets for the 2D and 3D cell approaches respectively, in which the path effects only consist of anelastic attenuation. For both controlled datasets, the locations of events and sites are based on CS15.4-600. In this way, we reproduced the same paths as in CS15.4-600. All components of the total residuals ( L2L(�� ⃗ x e ) , S2S(�� ⃗ x s ) , P2P(�� ⃗ x e , �� ⃗ x s ), B 0 e and WS 0 e,s ) of the controlled datasets were generated randomly. In both controlled datasets, L2L(�� ⃗ x e ) and S2S(�� ⃗ x s ) were randomly drawn  N(0, 0.4 2 )), while B 0 e and WS 0 e,s were randomly drawn from N(0, 0.2 2 ). For the 2D approach, the study region was divided into 20 patches, each with a distinct attenuation rate δ i , which was randomly drawn from N(0, 0.001 2 ) (Fig. 14a). P2P(�� ⃗ x e , �� ⃗ x s ) was computed as the sum of attenuation of all cells from the closest point to the site (Eq. 6). For the 3D approach, P2P(�� ⃗ x e , �� ⃗ x s ) was computed from 3D cells. In the top 10 km, each 1 km-thick layer was divided into 20 patches, with each formed cell assigned a distinct attenuation rate δ i , which was also drawn from N(0, 0.001 2 ) (Fig. 15). Below 10 km, the entire layer was one single patch, whose δ i was Fig. 10 Results of δ with the 3D cell approach at selected depths. Each row denotes one reduced dataset drawn from the same normal distribution. The path effects were calculated as the cumulative attenuation along a straight line from the hypocenter to the site (Eq. 6). After applying the VCM to both controlled datasets, we found that the attenuation patterns can be accurately recovered in both cases (Figs. 14,15). For most cells, the differences between the output and input δ i are very close to zero. The only non-zero differences occurred at a few cells on the edge, which were only sampled by a few paths. This test shows that the cell approaches are able to capture the anelastic attenuation when it is dominant. In the Cyber-Shake simulations, for relatively long periods (e.g., 5 s), the anelastic attenuation itself is quite small, thus we infer that other effects such as source radiation pattern and scattering lead to inaccurate estimation. In most regions, we have no knowledge of the genuine anelastic attenuation, thus this finding serves as an important reminder when interpreting cellspecific attenuation results based on empirical datasets, especially at long periods. In other words, δ i should be considered as effective attenuation combining multiple non-modeled effects, instead of the genuine anelastic attenuation.
Another issue with the cell-specific attenuation approach is the single-line assumption of wave-propagation path. To study the sensitivity of the cell approaches, we performed another test with the 3D controlled dataset. After computing the input P2P(�� ⃗ x e , �� ⃗ x s ) from the hypocenter to a site, we randomly shifted the hypocenters by a small distance on the rupture plane N(0, σ 1 2 ), resulting in a subtle difference between the path used in computing the input path effects and the path used for regression in the VCM. Even with a very small shift (σ 1 = 1 km), significant deviations between the output and input δ are observed (Fig. 16). As a result, it appears that the cell approaches are highly sensitive to the assumption of seismic wave-propagation paths. As mentioned in Sect. 3, all the events in CS15.4 are above magnitude 6. The rupture plane of a magnitude 6 earthquake typically has a length of about 30 km and a width of about12 km. For such large rupture planes, the point source assumption used in the VCM is no longer appropriate. In the near future, we will continue to investigate the more accurate representation of wave-propagation paths in CyberShake simulations and better ways to decompose the systematic effects so that the genuine path effects can be captured. This issue has little impact on empirical regressions, as most events have small magnitudes and can be treated as point sources. An issue arises from empirical datasets for larger magnitude events however, as we do not know the correct answer from the simulations, of whether path effects are due to source, wave reflections/refractions, or anelastic attenuation; the empirical determination of path effects should therefore probably include increased epistemic uncertainty. Unlike the cell approaches, the midpoint approach only considers the aggregated attenuations from one event to one site. As a result, the strongest attenuations within the basins are recovered. However, the midpoint approach also smooths out small-scale features in attenuation patterns by aggregating them, which results in large correlation lengths of  (Fig. 9). Since R-INLA does not allow specifying the lower-or upper-bound of the correlation length, we investigated finding alternative ways to constrain correlation length of P2P(�� ⃗ x e , �� ⃗ x s ) during the VCM regression. One way to quantify the spatial correlation of path effects is to apply the semivariogram model to path terms computed from mixed-effects regression (Walling and Abrahamson 2012). The range of values in the bestfitting semivariogram function denotes the lag distance beyond which path terms are no longer correlated. We first computed the path terms with the mixed-effects regression for CS15.4-1000 and performed a regression on the semivariogram function with a spherical model (Fig. 17a). The best-fitting range value, 1.209, was used as the fixed correlation length of P2P(�� ⃗ x e , �� ⃗ x s ) in the VCM regression. The updated pattern still captured the strongest attenuation within the basins and much weaker attenuation outside the basins (Fig. 17b). More importantly, the weak attenuations were constrained at locations with data. This approach is promising for broader applications but requires additional investigation using recorded data.

Future work
In this work, the spatial correlation of source and site effects were put into decomposed residual components from an ergodic GMM, while Landwehr et al. (2016) have spatially correlated GMM coefficients. Having spatially correlated GMM coefficients is more advantageous that it is a more flexible framework to account for spatially correlated effects, which may better capture effects that are not included in source, site, and path effects (e.g., directivity). On the other hand, the Landwehr et al. (2016) approach is more computationally expensive. Moreover, it did not take genuine path effects into account as the Fig. 16 Maps of differences between input and output δ for the 3D controlled dataset with the 3D cell approach and shifted hypocenters coefficients that account for distance scaling only considers an average distance attenuation for each event in all directions. Therefore, systematically comparing the results of different components between the approaches used in Landwehr et al. (2016) and in this paper is highly recommended. Another crucial finding from this work is that the path effects are very sensitive to the assumption of the wave travel path. As a result, it is very challenging to compute path effects from large earthquakes from both empirical and simulation datasets as their rupture planes may extend to several hundred kilometers and one single path assumption is not appropriate. Although the rupture planes of large earthquakes can be very large, they also often have patches of elevated slips that are responsible for the majority of seismic energy radiation (i.e., subevents). However, these issues cannot be investigated with empirical datasets due to the paucity of large repeated events in a given region, and the unknowable detailed characteristics of the ruptures. Therefore, one possible solution would be to conduct future ground-motion simulations that include large and small earthquakes along the same rupture planes and study how path effects from the large earthquake relate to combinations of path effects from many sub-events or smaller earthquakes along the same rupture plane. Once a general relationship is established, approximating the path effects of future large earthquakes by combining path effects from smaller events (either empirical or simulated) on the same rupture plane might be possible.

Conclusions
This study applied a novel technique, the VCM, to a CyberShake simulation dataset. By definition, CyberShake is a non-ergodic model, and the various inputs used for the simulations are known, making it an ideal testbed for better understanding and testing of the capabilities of the VCM. We found that the obtained site effects are consistent with the simulation datasets. Furthermore, they are consistent in trends with empirical datasets, suggesting they could be extracted using the VCM and used to improve future non-ergodic GMMs for common site coverage. The VCM fails to capture the simulated source effects because the dense linear distribution of self-similar events in the CyberShake dataset violates the isotropic assumption required by the Gaussian Process. In this case, one may compute the point-estimate of the event term instead. For the path effects, both the 2D and 3D cell approaches failed to capture the anelastic attenuation from the 3D velocity model, mainly because other unmodeled effects like source radiation pattern and scattering lead to inaccurate estimation for longer periods. Moreover, CyberShake only has large magnitude events, thus the large rupture planes and surface wave generation lead to rather complicated wave propagation effects. It is inappropriate to assume that seismic waves propagate along a straight line between two points. A more sophisticated representation of wave propagation is necessary to recover the anelastic attenuation pattern. This finite-fault effect was one of the motivations for the development of large-scale simulations in the first place. The midpoint approach tends to better recover the expected anelastic attenuation patterns, especially when fixed correlation length is carefully selected, as it does not make any assumption on the wave-propagation path. This study provides important guidance for future applications of VCM to CyberShake datasets and other simulation datasets. Due to the large discrepancy in magnitude range and earthquake density between CyberShake and empirical datasets, the specific issues encountered in the computation of source and path effects are unlikely to appear for empirical datasets. However, we expect that the general lessons learned will encourage further verifications from modelers, especially as the VCM is used to extrapolate ground motions to events much larger in magnitude than those empirically available. Last but not least, lessons learnt here will help further development of the technique to capture the genuine path effects for non-ergodic modelling development, both from simulations and as part of GMMs.