The Dabaghi and Der Kiureghian stochastic near-fault ground motion model requires information about the source, site, and source-to-site geometry, including directivity parameters. Directivity parameters entail often unavailable knowledge of the rupture geometry and hypocenter location. This article presents methods to randomize the directivity parameters required to simulate near-fault ground motions. A first procedure is proposed where only the contributing fault, earthquake magnitude, and site location are known. Possible rupture directivity conditions are accounted for by randomizing the rupture geometry and hypocenter location. For this purpose, new predictive models of the rupture geometry parameters are developed for shallow crustal earthquakes with magnitudes between 5.2 and 7.9. To allow validation of synthetic motions with NGA-West2 models, a second procedure randomizes the rupture geometry and both hypocenter and site locations. Results show a general agreement between the two methods.
Near-fault ground motions tend to be more complex and more variable than far-field ground motions and may possess distinct characteristics, such as the rupture directivity, hanging-wall (HW) effects, and the fling step, which may have profound influences on structural response. The rupture-directivity effect depends on the direction, length, and velocity of rupture propagation relative to a considered site (Somerville et al., 1997). When the rupture front propagates toward the site at a velocity nearly equal to the shear-wave velocity of the ground, all the seismic energy radiated from the fault rupture tends to arrive at the site in a single, large-amplitude, short-duration pulse—the so-called forward directivity pulse—in the fault-normal direction. If the rupture propagates away from the site, the site is in the backward directivity region and typically records a small amplitude, long-duration ground motion in the fault-normal direction.
The assessment of seismic risk and performance of structures located near active faults should account for the distinct near-fault characteristics of potential earthquakes. Given the scarcity of recorded near-fault ground motions from larger magnitude events and exhibiting near-fault effects, synthetic ground motions can be useful in seismic assessment studies in addition to or in place of recorded motions. Obviously, these synthetic ground motions must capture the important characteristics of near-fault ground motions and properly represent their natural variability for given earthquake source and site characteristics.
Several ground motion simulation techniques have been proposed. Deterministic physics-based (or source-based) techniques explicitly model the rupture process at the seismic source, propagate the resulting seismic waves, and determine the site response (Hartzell et al., 2005). Using a three-dimensional (3D) velocity model, these methods capture the local wave-propagation and basin effects. However, the generated synthetic motions typically lack high-frequency components beyond approximately 1 Hz. To address this shortcoming, high-frequency stochastic components are added to the physics-based synthetics (Graves and Pitarka, 2016). Physics-based approaches require a precise definition of the fault and rupture geometry and the shear-wave velocity of the medium, which often is not available to engineers.
Dabaghi and Der Kiureghian (2014, 2017, 2018) proposed a site-based stochastic model and simulation procedure to generate orthogonal horizontal components of synthetic near-fault ground motions for specified earthquake source and site characteristics. This method accounts for rupture directivity and produces pulse-like and non-pulse-like motions in accordance with their observed proportions among recorded motions. The model requires specification of the source (type-of-faulting , earthquake magnitude , depth to the top of the rupture plane ), the shear-wave velocity of the top 30 m of soil at the site (), and the source-to-site geometry (closest distance , directivity parameters and ). As illustrated in Figure 1, depending on the nature of the fault, the directivity parameter describes the length or width of the portion of the rupture that propagates between the hypocenter and the site in the direction of the slip, and the directivity parameter describes the angle in a horizontal plane between the fault rupture plane and the direction between the epicenter and the site, or the angle in a vertical plane between the fault rupture plane and the direction between the hypocenter and the site. Parameters and are used for strike-slip faulting, while and are used for dip-slip faulting (Somerville et al., 1997). The directivity parameters influence the probability that a ground motion is pulse-like and, if so, the amplitude and period of the pulse.
Ground-motion prediction equations (GMPEs) developed as part of the next-generation attenuation (NGA)-West2 project (Bozorgnia et al., 2014) do not explicitly account for the directivity effect, except those proposed by Chiou and Youngs (2014). More recently, five directivity models were developed based on the NGA-West2 ground motion database and numerical simulations of large earthquakes (Spudich et al., 2013, 2014). Comparisons of these models showed that their predictions have high variability and strong dependence on individual assumptions, especially for non-strike-slip faults (Spudich et al., 2013, 2014). These models are not considered sufficiently developed to be incorporated into the GMPEs and, therefore, further studies have been recommended (Spudich et al., 2013, 2014).
The near-fault simulation procedure of Dabaghi and Der Kiureghian (2018) requires as input the parameter set (, , , , , , ). To perform seismic design or assessment at a site with known location and , structural engineers can obtain information about the magnitude , the distance , and the source faults of earthquakes that contribute most to the hazard from deaggregation of probabilistic seismic hazard analysis, for example, as provided by the US Geological Survey (USGS, 2019) unified hazard tool. However, structural engineers are not likely to have information about the parameters , or even , which require knowledge of the rupture geometry as well as the location of the hypocenter.
Various models have been proposed to describe the scaling of rupture geometry parameters (e.g. length, width, aspect ratio, area) with magnitude and vice versa (Ellsworth, 2003; Goda et al., 2016; Hanks and Bakun, 2002; Leonard, 2010; Shaw, 2009; Thingbaijam et al., 2017; Wells and Coppersmith, 1994). These models employ various functional forms and are developed by regression analysis using diverse historical datasets of earthquakes from different tectonic regions. For a given , they can predict different rupture dimensions. With the exception of the study by Goda et al. (2016), these models predict the rupture length and width independently, without accounting for possible correlation between them. Moreover, they do not account for possible differences in the scaling of rupture geometry parameters between buried and surface ruptures.
Fewer studies are available for predicting . Gupta (2006) used the deterministic value , independent of ; Campbell et al. (2009) used a median estimate of that depends on , while Kaklamanos et al. (2011) estimated based on hypocentral depth, down-dip rupture width , and dip angle , assuming a deterministic hypocenter location at 60% down . Chiou and Youngs (2014) proposed a model that relates the mean of to and the type of faulting, and used it in their own NGA West2 GMPE. However, their model does not account for the variability of .
The objective of this article is to develop procedures to simulate the parameters , , , and , when they are unknown, to generate pairs of horizontal components of synthetic near-fault ground motions using the model and simulation procedure of Dabaghi and Der Kiureghian (2018). Two procedures are developed as described below.
The first procedure (P1) generates the parameters (, , , ) randomly for a specified site, given the type of faulting (), the magnitude () of the earthquake, the geometry of the contributing fault (location, strike, dip, maximum dimensions), and the location of the site. This is achieved by randomizing the fault rupture geometry, its location within the fault plane, and the hypocenter location within the rupture plane, based on their respective probability distributions. First, the rupture geometry is simulated using distributions of the rupture length , rupture width , and . These distributions are developed herein using regression analysis on a subset of earthquakes in the NGA-West2 database. Next, the location of the rupture plane is randomized within the fault plane. The hypocenter location along the strike and down-dip is then simulated based on the model developed by Mai et al. (2005). Finally, the coordinates of the site and of the simulated rupture plane and hypocenter are used to calculate and the directivity parameters and . This procedure is useful for structural engineers, who wish to simulate sets of acceleration time series for a particular site of interest for seismic design or assessment studies.
The second procedure (P2) randomly generates and for given , , , and and the dip angle of the fault. This procedure simulates the conditions present in the NGA-West2 dataset and is used to compare the synthetic ground motions with the GMPEs. This is achieved by randomizing the fault rupture geometry, the hypocenter location within the rupture plane, and the site location. The rupture geometry and hypocenter location are simulated as described above, while the site location is sampled from a uniform distribution along the locus of points at a distance from the simulated fault rupture.
After a brief review of the Dabaghi and Der Kiureghian (2014, 2017, 2018) stochastic near-fault ground motion model and simulation procedure, we develop predictive models for the rupture geometry parameters , , and in terms of the earthquake magnitude and type-of-faulting, and compare them to existing models. Our models are limited to shallow crustal earthquakes in active tectonic regions having magnitudes between 5.2 and 7.9. The hypocenter location model developed by Mai et al. (2005) and adopted in this article is then briefly reviewed. Next, procedure P1 is described and illustrated by simulating ground motions for a specified site given and and the geometry of the contributing fault. Procedure P2, which generates and for a random site for given , , , and and the dip angle of the fault, is then described. Finally, near-fault ground motions simulated using P2 are compared with the NGA-West2 GMPEs.
Stochastic near-fault ground motion model and simulation method
Dabaghi and Der Kiureghian (2014, 2017) developed a site-based parameterized stochastic model of near-fault ground motion in two orthogonal horizontal directions. The model accounts for the rupture directivity effect and considers both pulse-like and non-pulse-like motions. The model employs a modulated and filtered white-noise (MFW) process with time-varying filter parameters (Rezaeian and Der Kiureghian, 2008), and a modified version of the Mavroeidis and Papageorgiou (2003) pulse model. The model is defined in terms of physically meaningful parameters that represent the amplitude, duration, frequency content, and pulse characteristics of near-fault ground motions. The model was fitted to recordings taken from the Pacific Earthquake Engineering Research Center (PEER) NGA-West2 database to obtain empirical “observations” of the model parameters. The dataset used was limited to near-fault recordings from shallow crustal strike-slip, reverse, or reverse-oblique earthquakes in active tectonic regions with 5.5 ≤ ≤ 8.0 and recorded at sites located at < 31 km. Empirical predictive equations were constructed for the model parameters in terms of earthquake source and site characteristics, namely , , , , , , and . The variability of the model parameters and statistical correlations between them were estimated. Dabaghi and Der Kiureghian (2014, 2018) extended their model to generate suites of horizontal pairs of synthetic near-fault ground motions for specified earthquake scenarios. Figure A1 in Supplemental Appendix A shows a flowchart of this simulation procedure, where the empirical model of the probability of occurrence of the directivity pulse by Shahi and Baker (2014) is used to generate pulse-like and non-pulse-like ground motions in appropriate proportions. The resulting motions were shown to have realistic time series and to reproduce important features of recorded near-fault ground motions, including the directivity effect and natural variability. See the study by Dabaghi and Der Kiureghian (2014, 2018) for more details.
Rupture geometry
We wish to generate random realizations of the parameters (, , , ) according to probability distributions that are consistent with the physics of the earthquake rupture process and the available data. We start by developing new predictive models for , and with as the predictor variables, but consideration is also given to the dip angle of the fault and the seismogenic depth (thickness of the crust), if they are available. Differences in the scaling of rupture dimensions for buried () and surface () ruptures are also considered. The study is limited to crustal earthquakes in active tectonic regions.
Database of earthquake ruptures
A subset of records from the PEER NGA-West2 database (PEER, 2015) was used to fit the Dabaghi and Der Kiureghian (2014, 2018) model. The NGA-West2 database, which is limited to shallow crustal earthquakes in active tectonic regions, also provides the rupture geometry and dimensions for each earthquake. This database is used again in this study to fit the rupture geometry models. The selected subsets are those earthquakes whose rupture geometry and dimensions are determined by inversion of seismic waveforms. Included are 77 earthquakes with from strike-slip (SS), reverse (REV), and normal (N) faulting mechanisms (see Table B1 in Supplemental Appendix B for more details). Figure 2 is a scatter plot showing the distribution of earthquakes in the considered database with respect to moment magnitude and type-of-faulting. We note that data from normal and normal–oblique faults as well as data from reverse and reverse–oblique faults are combined throughout the analysis.
Distribution of earthquakes in the considered database with respect to moment magnitude and type of faulting.
Depth to top of rupture,
Data indicate that is zero-inflated, that is, the zero value is frequently observed, as can be seen in Figure 3. Thus, the value of is used to classify ruptures as surface , that is, the rupture extends to the Earth’s surface, or buried ), that is, the rupture does not reach the surface. We note that other references define a surface rupture differently. For example, Kagawa et al. (2004) define a surface rupture as having a clear surface dislocation and significant slip, and as a buried rupture otherwise.
Earthquake magnitude versus depth to the top of rupture in the considered database.
We employ a two-part model to represent the zero and non-zero values of . The first is a logistic regression model that predicts the probability that . The second is a linear regression model that predicts the mean and variability of the non-zero data. While developing these models, it was found that the data from normal faulting are not statistically different from the data obtained for strike-slip or reverse faulting. This could be due to the limited number of earthquakes with normal faults in the considered database. For this reason and following the study by Chiou and Youngs (2014), in this study, the data from normal faulting are combined with the data from strike-slip faulting to fit the models. The fitted logistic regression model is given by:
and is plotted in Figure 4 alongside the data used for fitting.
Logistic regression models for against used data.
To develop predictive equations for when , and similarly for and in the next section, the data are first transformed to ensure normality of the regression residuals. In all cases, the logarithmic transformation is used. The functional form of the predictive equation of a transformed model variable is initially constructed as:
where denotes the regression residual, assumed to be a normally distributed random variable with a zero mean and standard deviation , and , …, are the model parameters. , 1, and 2 correspond to SS, REV, and N faults, respectively. Step-wise linear regression is utilized to investigate the importance of the predictors, and the models are refitted accordingly. The models are evaluated based on the resulting value, residual plots, plots of fitted versus observed values, and normality of the residuals using the Shapiro–Wilk normality test (Shapiro and Wilk, 1965). A significance level of 5% is used for all statistical tests unless stated otherwise.
As previously mentioned, following the study by Chiou and Youngs (2014), data from N faults are combined with SS faults when developing the models. Analysis of the data for leads to the model:
where is measured in kilometers, and and are the expected value and standard deviation from the regression model, respectively. Moreover, to prevent the model from giving unreasonable predictions at lower magnitudes and in the tail, upper limits are imposed on . Based on observed maximum values in the database, the following upper bounds are considered:
Figure 5 shows the predictive model for given by Equation 3 and the upper bounds given by Equation 4, plotted against the data used for fitting. It can be seen in Figure 4 that the probability of having a buried rupture decreases with increasing magnitude and is larger for reverse faults than for strike-slip or normal faults. For buried ruptures in Figure 5, tends to decrease with magnitude, and reverse faults tend to have larger values than strike-slip and normal faults. The observed differences between fault types could be explained by the fact that strike-slip and normal faults tend to be steeper than reverse faults; thus, the rupture width of a reverse fault can extend further along the dip before rupturing the surface (or reaching the seismogenic depth). These observations are consistent with previous models, such as the one proposed by Chiou and Youngs (2014). The reader is cautioned against extrapolation of these models beyond the range of the data ().
Predictive models for with upper limits, against used data.
Rupture dimensions and
As mentioned before, calculating , , and for a given site requires knowledge of the geometry and coordinates of the fault rupture and the location of the hypocenter. The fault rupture is idealized as a rectangle with the aspect ratio varying from earthquake to earthquake. In this section, predictive models are developed for and . Similar to the study by Wells and Coppersmith (1994) and Thingbaijam et al. (2017), we differentiate between strike-slip, reverse and normal faulting because of the differences in their rupture mechanisms and dip angles. Possible influence on the rupture dimensions depending on the buried () or surface () nature of the rupture is also considered.
A procedure similar to that used for developing the linear regression model of is used. For the rupture length, no statistical difference is observed between surface and buried ruptures, so a predictive model that does not make this distinction is developed. The result is as follows:
where is in kilometers. Note that the rate of change with magnitude is not found to be affected by the type of faulting; only the intercept is dependent on the type of faulting. Figure 6 shows the fitted predictive model for together with the employed data. The model indicates that increases with , and that SS ruptures tend to be longer than REV and N ruptures. These are consistent with previous findings (Schwartz, 2018). Figure 6 also shows the widely used model for subsurface rupture length developed by Wells and Coppersmith (1994), referred to as WC94. The comparison shows that, in general, the developed models are similar to the WC94 models and are a good fit to the available data for .
Comparison of developed predictive models for with WC94 and employed data.
For the rupture width, a significant statistical difference is observed between the surface and buried rupture cases. Thus, the following models are developed:
where is in kilometers. Note again that the rate of change with magnitude is not statistically found to be affected by the type of faulting; only the intercept is dependent on the type of faulting. Figure 7 shows the fitted predictive models for , together with the data used to fit them. The models indicate that increases with , and that REV ruptures tend to be wider than SS ruptures. As for N faults, their scaling is similar to that of REV faults for buried ruptures, and close to (slightly above) that of SS faults for surface ruptures. The models also indicate that, when the rupture reaches the surface, the rate of change of with decreases, which is indicative of possible saturation due to the finite width of the seismogenic layer, while the rate of change of with is not affected (see Equations 5 to 7). This is consistent with the observation of Dalguer et al. (2008) that surface ruptures tend to have a larger aspect-ratio than buried ruptures. Again here, the larger width of REV surface ruptures compared to SS and N surface ruptures could be explained by the smaller dip angles of REV faults compared to the steeper SS and N faults. Figure 7 also shows the WC94 rupture width models, which distinguish between the three faulting mechanisms but not between surface and buried ruptures.
Comparison of developed predictive models for with WC94 and employed data.
Upper bounds on and are imposed based on seismological considerations to avoid unrealistic values. If the length of the fault is known, the upper bound is imposed. Otherwise is assumed, based on the largest observed rupture length (Schwartz, 2018). The upper bound depends on the width of the fault. If is not known, is assumed, where is the seismogenic depth, is the depth to the top edge of the fault, and is the dip of the fault in degrees. If or are not available, default values and km (Watts and Burov, 2003) are assumed. If the dip angle is not available, mean values based on faulting mechanism may be used: 85° for SS faulting, 40° for REV faulting, and 53° for N faulting. These mean values are obtained based on fault dip data in the considered database of shallow crustal earthquakes and are consistent with the values reported by Thingbaijam et al. (2017). Thus, for a given value of , the upper bound on the rupture width is calculated as:
Next, correlations between the log-transformed rupture geometry parameters (, ) are examined. These are estimated as the correlations between the residuals of the regression models (see Equation 2) and are listed in Table 1 along with their 95% confidence intervals (in brackets; Witte and Witte, 2016). The data indicate that there is a mild positive correlation of 0.244 (at 8% significance level) between the logarithms of the rupture length and width only for buried ruptures , while for surface ruptures, length and width are effectively uncorrelated. Note that when surface and buried ruptures are combined, the correlation between the logarithms of the rupture length and width is found to be 0.150, consistent with the correlation reported in the study by Goda et al. (2016). Further investigation shows that the correlation between and tends to decrease with magnitude, which explains why buried ruptures, which are usually associated with smaller magnitude events, show a positive correlation, while surface ruptures, which are usually associated with larger magnitude events, show close to zero correlation. Moreover, for buried ruptures, appears to be uncorrelated with both and . Later, for simulation purposes, we assume a correlation of 0.244 between and for buried ruptures and zero correlation among , , and for all other cases.
Estimated correlation matrix of regression residuals and their 95% confidence intervals (in brackets)
For
For
1
0.026 [−0.36, 0.41]
1
−0.028 [−0.30, 0.25]
−0.121 [−0.38, 0.16]
0.026 [−0.36, 0.41]
1
−0.028 [−0.30, 0.25]
1
0.244 [−0.03, 0.49]
−0.121 [−0.38, 0.16]
0.244 [−0.03, 0.49]
1
Simulation of rupture plane within fault plane
Given the coordinates of the top points of the fault, the fault length and strike are determined. In addition knowing the fault width and dip angle , the rectangular fault plane is defined. For any set of the simulated parameters (), the location of the rectangular rupture at depth within the fault plane is randomized assuming a uniform distribution within the fault plane.
Simulation of hypocenter location within rupture plane
The location of the hypocenter within the rupture plane has an important influence on the near-fault ground motion since it influences the directivity effect (Somerville et al., 1997). The location of the hypocenter is a required input for physics-based simulations (Graves and Pitarka, 2016; Hartzell et al., 2005; Pitarka et al., 2017) and it is also required to determine the directivity parameters and of the Dabaghi and Der Kiureghian (2018) model. To properly account for the variability of the generated near-fault ground motions, the simulation procedure should account for the variability in the hypocenter location. Some previous simulation studies, for example, Graves and Pitarka (2016) and Pitarka et al. (2017), have neglected the variabilities in the fault rupture area and hypocenter location.
Most previous studies assume the location of the hypocenter is uniformly distributed along the fault strike. The down-dip location is assumed to be either deterministic or to also follow a uniform distribution. For example, early studies in the CyberShake project placed hypocenters every 20 km along the strike (a discrete uniform distribution) and at a down-dip distance of 0.75 of the rupture width (Graves et al., 2011), while more recent studies placed them every 4.5 km along the strike and down-dip (Southern California Earthquake Center (SCEC), 2018). However, Mai et al. (2005) found that data from more than 50 earthquakes do not support the assumption of uniform distribution of the hypocenter location over the rupture plane. Their data indicated that the hypocenter is more likely to be near the center of the rupture plane than near the edges. They modeled the location of the hypocenter along the strike as a random variable following a truncated normal distribution with a mean of 0.5LR and a standard deviation of 0.23LR, while the down-dip location normalized by the rupture width was modeled using a truncated Weibull distribution. The estimated shape and scale parameters of the fitted Weibull distributions are, respectively, 0.626 and 3.921 for SS faults and 0.692 and 3.394 for REV and N faults. We adopt these models in this article, while assuming that the along-strike and down-dip locations of the hypocenter within the rupture plane are statistically independent.
Procedure P1: simulation for specified site and random rupture geometry and hypocenter location
We first consider the case where , , the fault geometry (coordinates of the top points, fault dimensions and , depth to the top edge of the fault , and dip angle ) and the site location are known, that is, procedure P1. In addition, the site is assumed to be known. As mentioned earlier, if some of the fault geometry parameters are unavailable, the default values , km, km and for SS faulting, 40° for REV faulting, and 53° for N faulting are used. The following procedure is employed to simulate , and the directivity parameters and :
Given and , compute the probability of obtaining from Equation 1. Randomly generate the number 0 or 1 with probabilities and , respectively. If number 0 is generated, , otherwise, .
If , that is, the case of surface rupture, simulate and as normally distributed uncorrelated random variables with means and standard deviations computed from Equations 5 and 6, respectively. If , that is, the case of buried rupture, simulate , , and as jointly normally distributed random variables with means and standard deviations estimated from Equations 3, 5, and 7, respectively, and a correlation coefficient of 0.244 between and . If any of the simulated variables exceeds its specified upper bound, repeat step 2 until all upper bounds are satisfied.
Given the geometry of the fault plane and the simulated rupture parameters (, , and ), uniformly sample the location of the rupture plane within the fault plane and calculate the corresponding coordinates of the edges of the rupture plane.
Given , and , simulate the location of the hypocenter along the strike and down-dip according to the Mai et al. (2005) model. Use the coordinates of the rupture plane and to calculate the coordinates of the hypocenter.
Use the coordinates of the site and the simulated coordinates of the rupture plane and the hypocenter to determine , , and .
Each simulation results in a different set of , , , and rupture and hypocenter locations, thus in different values of , , and . This process is repeated to generate any desired number of rupture directivity scenarios for a specified set of input values , , fault geometry, and site location). In addition, given of the site, the resulting (, , , , , , ) scenarios are provided as input to the near-fault simulation procedure of Dabaghi and Der Kiureghian (2018) to generate the desired number of synthetic near-fault ground motions. Figure 8 presents simplified flowcharts of the Dabaghi and Der Kiureghian (2018) simulation procedure (left) and of simulation procedure P1 (center), and the relation between them. A more detailed flowchart of P1 can be found in Figure A2 of Supplemental Appendix A.
Simplified flowcharts of the Dabaghi and Der Kiureghian (2014, 2018) near-fault (NF) ground motion (GM) simulation procedure (left), simulation procedure P1 (center), and simulation procedure P2 (right). More detailed flowcharts of each of these simulation procedures can be found in Supplemental Figures A1 to A3.
Example application
The proposed simulation procedure is used to generate pairs of horizontal components of near-fault ground motions for an earthquake event at a site selected from the CyberShake platform (SCEC, 2018) in downtown Los Angeles (LADT) having latitude 34.05204, longitude −118.25713, and m/s. The 2008 edition of the USGS hazard model together with its unified hazard tool is used to determine the earthquakes and faults that contribute most to the hazard at the selected LADT site. For example, at a spectral period of 2 s and for a return period of 2475 years, the deaggregation results show a major contribution from earthquakes having large magnitudes () and occurring at short distances (< 10 km) from the site. One of the largest contributors to the hazard is an earthquake with on the Upper Elysian Park reverse fault. This earthquake scenario is used to illustrate our proposed simulation procedure. Information about the geometry of the Upper Elysian Park fault is extracted from the CyberShake platform, which is based on Version 2.0 of the Uniform California Earthquake Rupture Forecast, UCERF2.0 (Field et al., 2009). The extracted information includes the coordinates of the endpoints of the top edge of the fault, from which the fault length km is calculated as the horizontal distance between the two points, that is, by idealizing the fault as a rectangular plane. The above reference also provides km, km, and .
Given information about the fault type (REV) and geometry, and the site location, procedure P1 is used to generate, for , 100 different realizations of , , , and rupture and hypocenter locations within the fault plane, which are in turn used to calculate the corresponding , , and values. Table 2 summarizes the resulting ranges of some of these simulated variables and Figure 9 illustrates selected realizations of the rupture geometry and hypocenter locations.
Ranges of simulated rupture geometry parameters and calculated distance and directivity parameters
Simulated parameter
(km)
(km)
(km)
(km)
(km)
(°)
Min
7.9
5.2
3.1
6.4
1.7
0.4
Max
18.4
20.3
9.3
10.8
17.0
16.8
Selected realizations of rupture geometry (gray rectangular surface) and hypocenter location (asterisk) along strike and down-dip for an earthquake on the Upper Elysian Park fault (dashed rectangle).
For the 100 different sets of simulated rupture realizations, the corresponding parameters , , , and together with (REV), , and m/s are used as input to the simulation procedure of Dabaghi and Der Kiureghian (2018) to obtain 100 pairs of horizontal components of near-fault ground motions. The 5% damped RotD50 pseudo-acceleration response spectra (Boore, 2010) of the simulated motions and their median and median plus and minus one logarithmic standard deviations are illustrated in Figure 10a. These simulations, which cover a range of directivity configurations (see Table 2), consist of 38 pulse-like and 62 non-pulse-like motions. For the former set, the pulse period ranges from 0.44 to 2.44 s. Over this period range, the pulse-like motions tend to have larger amplitudes than non-pulse-like motions. Moreover, Figure 10a highlights one pulse-like simulation (thick dark gray line) with s and resulting from a forward directivity scenario ( km), and one non-pulse-like simulation (thick light gray line) resulting from a backward directivity scenario ( km). Note that the spectral amplitudes of the pulse-like motion are amplified near the period of the pulse. The acceleration, velocity, and displacement time series of the two highlighted synthetic ground motions are shown in Figure 10b.
(a) RotD50 spectra of 100 simulated motions (38 pulse-like and 62 non-pulse-like) at the LADT site ( m/s) due to an earthquake with occurring on the Upper Elysian Park fault; median and median plus and minus one logarithmic standard deviations spectra, and examples of one pulse-like motion (from a scenario with km and with s) and one non-pulse-like motion (from a scenario with km) are shown. (b) Acceleration, velocity, and displacement time series of the two horizontal components of the pulse-like (dark gray) and non-pulse-like (light gray) ground motions that are highlighted in (a).
Procedure P2: simulation for random site, rupture geometry and hypocenter location
To compare the simulation results with NGA-West2 GMPEs, we consider the case where, in addition to , , and the of the site, and are also known, that is, procedure P2. Unless available, the default values are used for the fault geometry parameters , , , and , similar to P1. The following procedure is developed to simulate rupture directivity conditions (, ):
Given , , and , simulate and as normally distributed uncorrelated random variables with means and standard deviations estimated from Equations 5 and 6, respectively, if , or as jointly normal random variables with means and standard deviations estimated from Equations 5 and 7, respectively, and correlation coefficient 0.244, if . Follow a procedure similar to step 2 in P1 to ensure sampled and values are within the specified upper bounds.
Given , , and , simulate the location of the hypocenter within the rupture plane similar to step 4 in P1. Use and to calculate the coordinates of the hypocenter.
Given the simulated rupture geometry, we determine a “racetrack” at a distance from the simulated fault rupture using the method developed by Brian Chiou (2018, personal communication). The “racetrack” represents the locus of all points on the ground surface that are located at a distance from the simulated rupture plane. Use the uniform distribution to sample a single site on the racetrack.
Use the simulated rupture geometry and hypocenter and site locations to determine the directivity parameters and .
Each simulation results in different , , and hypocenter and site locations and, therefore, different values of and . This process is repeated to generate any desired number of rupture directivity realizations for specified , , , , and . This allows a fair comparison of the synthetic motions with the NGA-West2 ground motion prediction equations. Figure 8 (right) presents a simplified flowchart of simulation procedure P2, and its relation with the Dabaghi and Der Kiureghian (2014, 2018) simulation procedure. A more detailed flowchart of P2 is presented in Figure A3 in Supplemental Appendix A.
Example application
The proposed simulation procedure P2 is used to generate 600 pairs of orthogonal horizontal components of near-fault ground motions for different hypothetical earthquake scenarios with selected , ,, , and values and random directivity conditions. For each simulation, , , and hypocenter and site locations are generated following the procedure described in the previous section. The resulting directivity variables are used as input to the simulation procedure of Dabaghi and Der Kiureghian (2018) to obtain the 600 pairs of orthogonal horizontal components of the ground motion.
The statistics of elastic response spectra of the generated synthetic ground motions are compared to those of the NGA-West2 GMPEs (Bozorgnia et al., 2014), and by that indirectly to recorded ground motions. We use a weighted average of the five NGA-West2 GMPEs developed by Abrahamson et al. (2014), Boore et al. (2014), Campbell and Bozorgnia (2014), Chiou and Youngs (2014), and Idriss (2014), respectively, referred to as ASK14, BSSA14, CB14, CY14, and I14. Following the recommendation by Rezaeian et al. (2014), all models are assigned a weight of 2/9 except I14, which is assigned a weight of 1/9. The weighted average of the five NGA-West2 GMPEs is denoted the “NGA-West2 model.”
The above GMPEs have limitations in the near-fault region. First, they do not require directivity parameters as input (except CY14) and, therefore, do not differentiate between forward and backward directivity sites. Thus, effectively, they represent random directivity conditions. Moreover, they were calibrated to broader ranges of magnitude and distance than the Dabaghi and Der Kiureghian (2018) model, which is limited to records from large magnitude earthquakes () at short distances ( km). Therefore, the NGA-West2 GMPEs are not specifically calibrated to near-fault scenarios.
The synthetic ground motions are compared with the GMPEs for several hypothetical earthquake scenarios occurring on vertical SS () or dipping REV () faults, having magnitudes Mw = 6.5, 7, or 7.5, and ZTOR = 0 or 3 km. The considered sites have Vs30 = 360, 525, or 760 m/s and are located at RRUP = 5, 10, or 20 km from the fault rupture. Simulations were conducted for a total of 34 scenarios that are reported in the study by Daoud et al. (forthcoming). For each of the considered scenarios, statistics of the 5% damped pseudo-acceleration response spectra of the 600 synthetic motions are compared with those described by the NGA-West2 model. The compared statistics are the median and median plus and minus one logarithmic standard deviation levels for the RotD50 horizontal component.
The NGA-West2 GMPEs require additional input parameters not required by the Dabaghi and Der Kiureghian model. For instance, some GMPEs need as input other distance parameters in addition to or instead of . These distance parameters include the Joyner–Boore distance (i.e. the closest distance to the surface projection of the coseismic rupture), the horizontal distance from the top of the rupture measured perpendicular to the fault strike, and the horizontal distance off the end of the rupture measured parallel to the strike. Some GMPEs also need as input the depth of the hypocenter or the width of the rupture. For each earthquake scenario, these additional input parameters are calculated for each realization of the rupture geometry and hypocenter and site locations, and their values are then used for comparison of the synthetic motions with the GMPEs.
Moreover, the NGA-West2 GMPEs (except I14) differentiate between sites located on the HW, that is, on the down-dip side of the top of rupture, and sites located on the footwall (FW). The HW effect refers to the increase in ground motion amplitudes at short periods observed at sites located at short distances on the HW side of the rupture, as compared to sites located at the same distance but on the FW side (Donahue and Abrahamson, 2014). This effect is accounted for in the NGA-West2 GMPEs using an input HW factor, herein referred to as (= 1 for HW sites; = 0 for FW sites). The simulation procedure of Dabaghi and Der Kiureghian does not explicitly distinguish between HW and FW sites. Thus, the simulated motions are compared with the GMPEs for random HW–FW configurations. Default values are assigned for the other input parameters that remain unspecified in the various GMPEs.
Figure 11 shows the comparison for two of the considered near-fault scenarios for random HW–FW configurations, one for an SS fault () and one for an REV fault (). In addition to the median and median plus/minus one standard deviation of the synthetic ground motion spectra, we show the median spectra for the subsets of pulse-like (P) and non-pulse-like (NP) synthetic motions. The generated synthetic ground motions show good agreement with the NGA-West2 GMPEs. At most periods and for both the strike-slip () and reverse () faulting mechanisms, the median spectrum of all the simulated motions falls within the range spanned by the median spectra of the five GMPEs. Moreover, the median level of the simulated motions falls within the median plus and minus one standard deviation levels predicted by the NGA-West2 model. For strike-slip faults, the synthetic motions tend to predict larger spectral ordinates at periods between 0.1 and 0.3 s and at longer periods (greater than 1 s) and to predict smaller spectral ordinates at other periods; see Figure 11 (left). However, the differences are not large. These observations are mostly consistent among all strike-slip scenarios. For the particular reverse faulting scenario shown in Figure 11 (right), the synthetic motions tend to predict larger spectral ordinates only at longer periods (greater than 1 s) and smaller spectral ordinates at all other periods. While these differences are larger than the ones observed for the strike-slip scenario, they are still not significant. Other reverse faulting scenarios exhibit slightly different trends, whereby some have more significant differences between the simulated motions and the GMPEs, while others show better agreement between the two. For events with , the simulations tend to predict spectral ordinates that are in overall good agreement with the NGA-West2 model at longer periods, but predict smaller spectral ordinates at shorter periods (less than 1 s).
Median and median plus and minus one logarithmic standard deviation of RotD50 spectra of 600 synthetic motions; median spectra of subsets of pulse-like (P) and non-pulse-like (NP) synthetic motions; median and median plus and minus one logarithmic standard deviation spectra predicted by a combination of five NGA-West2 GMPEs, and median spectra predicted by each of the five NGA-West2 GMPEs for two different earthquake scenarios.
It can be observed from Figure 11 that, for both fault types, the median spectrum of the subset of simulated pulse-like motions has considerably larger amplitudes than the median spectrum for all simulated motions as well as the median spectrum of the NGA-West2 model for periods greater than 1 s. This is due to the forward directivity effect, which is present in the latter spectra only in a weighted average sense. On the other hand, the median spectrum of the subset of simulated non-pulse motions has smaller amplitudes than the median spectrum of all simulated motions and smaller values than the median spectrum of the NGA-West2 model at periods 0.5–3 s. This is due to the absence of long-period pulses in these simulated motions. Similar trends are observed for the other scenarios.
Figure 12a to c show the effect on the median spectra of the simulated ground motions and the NGA-West2 model of varying , , and , respectively, while keeping other scenario parameters fixed. It is observed in Figure 12a that the effect of varying is not captured similarly by the synthetic ground motions and the NGA-West2 model. For the NGA-West2 model, as increases, spectral amplitudes increase at all periods, whereas for the synthetic motions spectral amplitudes increase with at longer periods, while magnitude saturation is observed at lower periods. On the contrary, Figure 12b shows that, as decreases, the spectral amplitudes increase at all periods for both the NGA-West2 model and the synthetic motions. Moreover, for the scenarios considered in Figure 12b, the period-dependent differences between the NGA-West2 model and the synthetic motions are similar to the differences observed in Figure 11. Finally, it is observed in Figure 12c that the effect of varying is not captured similarly by the NGA-West2 model and the synthetic ground motions. For the NGA-West2 model, as decreases spectral amplitudes increase at all periods, but more at the longer periods. For the synthetic motions, the effect of is negligible at shorter periods and spectral amplitudes increase as decreases only at periods longer than about 0.5 s. The trends that the synthetic motions show with at lower periods are comparable to the ones observed by CB14. Since both the simulated motions and the NGA-West2 model are based on fitting assumed models to data of recorded motions, it is not possible to ascertain as to which of the two produces more accurate trends. Nevertheless, it is noted that the differences between the two models and their trends are relatively small. And although some differences are observed in the median spectral shapes between simulated motions and the NGA-West2 model, similarities are observed between the spectral shapes of simulated motions and those of the ASK14 and BSSA14 GMPEs, especially at shorter periods.
Median of the RotD50 spectra for 600 synthetic motions versus median spectra predicted by the NGA-West2 model (a) for scenarios having the same , and values but different values, (b) for scenarios having the same , and values but different values, and (c) for scenarios having the same , and values but different values.
Conclusion
The site-based stochastic model and simulation procedure proposed by Dabaghi and Der Kiureghian (2014, 2017, 2018) requires information about the source, site, and source-to-site geometry, namely the parameters , , , , , , and . The resulting simulations account for the near-fault rupture directivity effect and capture the natural variability of real ground motions. They also make the crucial distinction between pulse-like and non-pulse-like ground motions. However, the directivity parameters and entail information about the rupture geometry and hypocenter location that may not be available to users of the procedure, particularly design engineers. In this article, a procedure to simulate the rupture geometry and hypocenter location and thus generate synthetic motions with random directivity effect for a given earthquake scenario is proposed.
Using data from a subset of shallow crustal earthquakes in active tectonic regions with in the NGA West2 database, new predictive models are developed for the rupture geometry parameters (, , and ) as a function of and . These models are found to be consistent with existing models. Important features of the developed models include: a two-part model consisting of a logistic regression model and a linear regression model to account for the zero-inflation in values; account for the variability and correlations between the rupture geometry parameters; and models for and that distinguish between buried ruptures () and surface ruptures (). These models reveal several important observations about the scaling of rupture geometry with magnitude. For example, the models indicate that the scaling of the rupture width with magnitude is different for buried and surface ruptures, while the scaling of the rupture length with magnitude is not affected by .
The first simulation procedure (P1) developed in this article considers the case where only the contributing fault, earthquake magnitude, and site location and characteristics are known. This procedure is appropriate for design or assessment situations. Using Monte Carlo simulation, random rupture geometries and hypocenter locations are generated according to their probability distributions and are used to calculate the corresponding randomized rupture directivity conditions. The location of the rupture within the fault plane is assumed to follow a uniform distribution, while the hypocenter location is assumed to follow the distributions proposed by Mai et al. (2005). The procedure is illustrated by simulating near-fault ground motions at a site located in downtown Los Angeles.
The second simulation procedure (P2) takes as input the type-of-faulting, the earthquake magnitude, the source-to-site distance, and the site characteristics. In this case, the rupture geometry as well as both the hypocenter and site locations are randomized. This procedure allows comparison with existing GMPE models. The procedure is illustrated for several earthquake scenarios defined by the set of parameters , , , , and . Results show general agreement between the two, but also some differences. For instance, synthetic ground motions are generally consistent with the GMPEs for random HW–FW configurations, but they are not able to capture the HW effect, though this effect is rather small. On the contrary, the synthetic motions explicitly account for the directivity effect, while GMPEs account for this effect in an average sense.
The proposed predictive models of rupture geometry and simulation procedure P1 are necessary for structural engineers who want to generate synthetic motions to perform seismic design or assessment studies at a near-fault site with known location and . They only need information about the magnitude () and the source fault of the earthquakes that contribute most to the hazard at their site of interest. This information can be obtained from deaggregation of probabilistic seismic hazard analysis results.
Supplemental Material
sj-pdf-1-eqs-10.1177_87552930211030941 – Supplemental material for Simulation of near-fault ground motions for randomized directivity parameters
Supplemental material, sj-pdf-1-eqs-10.1177_87552930211030941 for Simulation of near-fault ground motions for randomized directivity parameters by Yara Daoud, Mayssa Dabaghi and Armen Der Kiureghian in Earthquake Spectra
Footnotes
Acknowledgements
We thank Ibrahim Alameddine and Tom Shantz for their insightful feedback and suggestions, and Brian Chiou for providing his racetrack generation codes. We also thank the anonymous reviewers whose constructive comments improved the paper.
Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the State of California through the Transportation System Research Program of the Pacific Earthquake Engineering Research Center (PEER). Any opinions findings, and conclusion or recommendations expressed in this material are those of the author(s) and do not necessarily reflect those of the funding agency.
ORCID iD
Mayssa Dabaghi
Data and resources
The rupture geometry data used to develop the models in this study can be obtained from the NGA-West2 flatfiles ().
Supplemental material
Supplemental material for this article is available online.
References
1.
AbrahamsonNASilvaWJKamaiR (2014) Summary of the ASK14 ground motion relation for active crustal regions. Earthquake Spectra30(3): 1025–1055.
2.
BooreDM (2010) Orientation-independent, nongeometric-mean measures of seismic intensity from two horizontal components of motion. Bulletin of the Seismological Society of America100(4): 1830–1835.
3.
BooreD. M.StewartJ. P.SeyhanE.AtkinsonG. M., 2014. Nga-west2 equations for predicting pga, pgv, and 5% damped psa for shallow crustal earthquakes, Earthquake Spectra30(3), 1057–1085.
4.
BozorgniaYAbrahamsonNAAtikLAAnchetaTDAtkinsonGMBakerJWBaltayABooreDMCampbellKWChiouBSJDarraghRDaySMDonahueJLGravesRWGregorNHanksTCIdrissIMKamaiRKishidaTKottkeAMahinSARezaeianSRowshandelBSeyhanEShahiSKShantzTSilvaWSpudichPStewartJPWatson-LampreyJWooddellKYoungsR (2014) NGA-West2 research project. Earthquake Spectra30(3): 973–987.
5.
CampbellKWAbrahamsonNPowerMChiouBSJBozorgniaYShantzTRobleeC (2009) Next Generation Attenuation (NGA) project: Empirical ground motion prediction equations for active tectonic regions. In: Proceedings of the 6th international conference on urban earthquake engineering, Tokyo, Japan, 3–4 March.
6.
CampbellKWBozorgniaY (2014) NGA-West2 ground motion model for the average horizontal components of PGA, PGV, and 5% damped linear acceleration response spectra. Earthquake Spectra30(3): 1087–1115.
7.
ChiouBSJYoungsRR (2014) Update of the Chiou and Youngs NGA model for the average horizontal component of peak ground motion and response spectra. Earthquake Spectra30(3): 1117–1153.
8.
DabaghiMDer KiureghianA (2014) Stochastic modeling and simulation of near-fault ground motions for performance-based earthquake engineering. PEER report no. 2014/20, December. Berkeley, CA: University of California, Berkeley.
9.
DabaghiMDer KiureghianA (2017) Stochastic model for simulation of near-fault ground motions. Earthquake Engineering & Structural Dynamics46(6): 963–984.
10.
DabaghiMDer KiureghianA (2018) Simulation of orthogonal horizontal components of near-fault ground motion for specified earthquake source and site characteristics. Earthquake Engineering & Structural Dynamics47: 1369–1393.
11.
DalguerLAMiyakeHDaySMIrikuraK (2008) Surface rupturing and buried dynamic-rupture models calibrated with statistical observations of past earthquakes. Bulletin of the Seismological Society of America98(3): 1147–1161.
12.
DaoudYDabaghiMDer KiureghianA (forthcoming) Chapter 3: Simulation of near-fault ground motions for randomized directivity parameters. In: PEER Report: Improvements to the Stochastic Modeling and Simulation of Ground Motions for Use in PBEE.
EllsworthWL (2003). Appendix D: Magnitude and area data for strike slip earthquakes in: US Geological Survey Open File Report no. 03-214: Earthquake Probabilities in San Francisco Bay Region: 2002–2031, by Working Group on California Earthquake Probabilities, Menlo Park, CA.
15.
FieldEHDawsonTEFelzerKRFrankelADGuptaVJordanTHParsonsTPetersenMDSteinRSWeldonRWillsCJ (2009) Uniform California earthquake rupture forecast, version 2 (UCERF 2). Bulletin of the Seismological Society of America99(4): 2053–2107.
16.
GodaKYasudaTMoriNMaruyamaT (2016) New scaling relationships of earthquake source parameters for stochastic tsunami simulation. Coastal Engineering Journal58(3): 1650010-1–1650010-40.
17.
GravesRPitarkaA (2016) Kinematic ground-motion simulations on rough faults including effects of 3D stochastic velocity perturbations. Bulletin of the Seismological Society of America106(5): 2136–2153.
18.
GravesRJordanTHCallaghanSDeelmanEFieldEJuveGKesselmanCMaechlingPMehtaGMilnerKOkayaDSmallPVahiK (2011) CyberShake: A physics-based seismic hazard model for Southern California. Pure and Applied Geophysics168(3–4): 367–381.
19.
GuptaID (2006) Defining source-to-site distances for evaluation of design earthquake ground motion. In: Proceedings of the 13th symposium on earthquake engineering, Indian Institute of Technology Roorkee, Roorkee, India, 18–20 December.
20.
HanksTCBakunWH (2002) A bilinear source-scaling model for M-log a observations of continental earthquakes. Bulletin of the Seismological Society of America92(5): 1841–1846.
21.
HartzellSGuatteriMMaiPMLiuP-CFiskMR (2005) Calculation of broadband time histories of ground motion, Part II: Kinematic and dynamic modeling using theoretical Green’s functions and comparison with the 1994 northridge earthquake. Bulletin of the Seismological Society of America95(2): 614–645.
22.
IdrissIM (2014) An NGA-West2 empirical model for estimating the horizontal spectral values generated by shallow crustal earthquakes. Earthquake Spectra30(3): 1155–1177.
23.
KagawaTIrikuraKSomervillePG (2004) Differences in ground motion and fault rupture process between the surface and buried rupture earthquakes. Earth, Planets and Space56(1): 3–14.
24.
KaklamanosJBaiseLGBooreDM (2011) Estimating unknown input parameters when implementing the NGA ground-motion prediction equations in engineering practice. Earthquake Spectra27(4): 1219–1235.
25.
LeonardM (2010) Earthquake fault scaling: Self-consistent relating of rupture length, width, average displacement, and moment release. Bulletin of the Seismological Society of America100(5A): 1971–1988.
26.
MaiPMSpudichPBoatwrightJ (2005) Hypocenter locations in finite-source rupture models. Bulletin of the Seismological Society of America95(3): 965–980.
27.
MavroeidisGPPapageorgiouAS (2003) A mathematical representation of near-fault ground motions. Bulletin of the Seismological Society of America93(3): 1099–1131.
PitarkaAGravesRIrikuraKMiyakeHRodgersA (2017) Performance of Irikura recipe rupture model generator in earthquake ground motion simulations with Graves and Pitarka hybrid approach. Pure and Applied Geophysics174(9): 3537–3555.
30.
RezaeianSDer KiureghianA (2008) A stochastic ground motion model with separable temporal and spectral nonstationarities. Earthquake Engineering & Structural Dynamics37(13): 1565–1584.
31.
RezaeianSPetersenMDMoschettiMPPowersPHarmsenSCFrankelAD (2014) Implementation of NGA-West2 ground motion models in the 2014 U.S. National Seismic Hazard Maps. Earthquake Spectra30(3): 1319–1333.
32.
SchwartzDP (2018) Review: Past and future fault rupture lengths in seismic source characterization—The long and short of it. Bulletin of the Seismological Society of America108(5A): 2493–2520.
33.
ShahiSKBakerJW (2014) An efficient algorithm to identify strong-velocity pulses in multicomponent ground motions. Bulletin of the Seismological Society of America104(5): 2456–2466.
34.
ShapiroSSWilkMB (1965) An analysis of variance test for normality (complete samples)†. Biometrika52(3–4): 591–611.
35.
ShawBE (2009) Constant stress drop from small to great earthquakes in magnitude-area scaling. Bulletin of the Seismological Society of America99(2A): 871–875.
36.
SomervillePGSmithNFGravesRWAbrahamsonNA (1997) Modification of empirical strong ground motion attenuation relations to include the amplitude and duration effects of rupture directivity. Seismological Research Letters68(1): 199–222.
SpudichPBaylessJRBakerJWChiouBSJRowshandelBShahiSKSomervillePG (2013) Final report of the NGA-West2 directivity working group. PEER report no. 2013/09, May. Berkeley, CA: University of California, Berkeley.
39.
SpudichPRowshandelBShahiSKBakerJWChiouBSJ (2014) Comparison of NGA-West2 directivity models. Earthquake Spectra30(3): 1199–1221.
40.
ThingbaijamKKSMaiPMGodaK (2017) New empirical earthquake source-scaling laws. Bulletin of the Seismological Society of America107(5): 2225–2246.
WattsABBurovEB (2003) Lithospheric strength and its relationship to the elastic and seismogenic layer thickness. Earth and Planetary Science Letters213(1–2): 113–131.
43.
WellsDLCoppersmithKJ (1994) New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement. Bulletin of the Seismological Society of America84(4): 974–1002.
Please find the following supplemental material available below.
For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.
For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.