Reduction of Bias and Uncertainty in Regional Seismic Site Amplification Factors for Seismic Hazard and Risk Analysis

Abstract

Site amplification factors in National Building Codes are typically specified as a function of the average shear wave velocity over the first 30 m (V_s30) or site class (A, B, C, D and E) for defined ranges of V_s30 and/or ranges of depth to bedrock. However, a single set of amplification factors may not be representative of site conditions across the country, introducing a bias in seismic hazard and seismic risk analyses. This is exemplified by significant differences in geological settings between East and West coast locations in North America. Western sites are typically characterized by lower impedance contrasts between recent surface deposits and bedrock in comparison to Eastern sites. In North America, site amplification factors have been derived from a combination of field data on ground motions recorded during West Coast earthquakes and numerical models of site responses that are meant to be representative of a wide variety of soil profiles and ground motions. The bias on amplifications and their impact on seismic hazards is investigated for the Montreal area, which ranks second for seismic risks in Canada in terms of population and hazard (PGA of 0.25 g for a 2475 years return period). Representative soil profiles at several locations in Montreal are analyzed with 1-D site response models for natural and synthetic ground motions scaled between 0.1 to 0.5 g. Since bedrock depths are typically shallow (<30 m) across the island, bedrock shear wave velocities have a significant influence on the impedance contrast and amplifications. Bedrock shear wave velocity is usually very variable due to the differences in rock formations, level of weathering and fracturing. The level of this uncertainty is shown to be greatly decreased when rock quality designation (RQD) data, common information when bore hole data are logged, is available since it is highly correlated with both shear and compression wave velocities. The results are used to derive region-specific site amplification factors as a function of both V_s30 and site fundamental frequency and compared to those of the National Building Code of Canada (2015). The results of the study indicate that there are large uncertainties associated with these parameters due to variability in soil profiles, soil properties and input seismic ground motions. Average and confidence intervals for the mean and for predictions of amplification factors are calculated for each site class to quantify this uncertainty. Amplifications normalized relative to class C are obtained by accounting for the correlation between site class amplifications for given ground motions. Non-linearity in the analysis of equivalent linear 1-D site response is taken into account by introducing the non-linear G/Gmax and damping ratios curves. In this method, it is assumed that the shear strain compatible shear modulus and damping ratio values remains constant throughout the duration of the seismic excitation. This assumption is not fully applicable to a case when loose saturated soil profile undergo heavy shaking (PGA > 0.3 g). In this study, all simulations with input motion PGA >0.3 g have been performed by using the EL method instead of the NL method considering that cohesive soils (clay and silt) at Montreal sites are stiff and cohesionless soils (sand and gravel) are considerably dense. In addition, the field and laboratory data required to perform NL analyses are not currently available and may be investigated in future works.

1. Introduction

The region of Montreal, Canada, is the second city at risk in Canada for its population and seismic hazard [1]. However, no damaging events have been recorded to provide insight on site response and structural performance. The most significant recent event in Montreal occurred in 1988 (M 5.7) in the Saguenay region, approximately 350 km NE of the city and caused structural damaged to the Montreal East city hall built on a 17 m soft clay layer that may have amplified the ground motions [2,3]. Zones with similar soft soil are found in the downtown area and along the south east shore of the island [4,5]. Site effects that are currently used in design codes may not reflect the amplification to be expected for site conditions and ground motions that are typical to eastern Canada as well as the uncertainty associated with these factors. Given the lack of historical data, an approach based on simulation of ground motions is proposed here to investigate site effects.

The equivalent linear 1-D ground response analysis method (later named 1-D EL) is often used to estimate site effects and provides the seismic motion at the ground surface of the soil column. For this purpose, the site amplification factor for a given period is defined as the ratio of the spectral acceleration computed at the ground surface to the spectral acceleration of the input motion applied at the rock interface. The site amplification factors of the 2015 National Building Code of Canada (NBCC-2015) are specified for site classes A, B, C, D and E (as defined in NHERP, 1994) as a function of structural periods. They are largely derived from analyses performed by Choi and Stewart [6] using 1828 recordings from 154 shallow crustal earthquakes and provide factors for NEHRP V_s30 soil classes inferred from geotechnical data at 209 strong motion stations. They presented the 1-D dynamic response analyses for typical soil profiles for rock motions scaled to 0.16 g using the equivalent linear model of SHAKE2000 [7] and developed an empirical relation between amplification of rock motions and site periods. In our study, results from the most recent seismic surveys are used to characterize shear wave velocity profiles in Montreal for soil classes B, C and D and to estimate site amplification factors from the 1-D ground response analysis using Strata [8] for a range of ground motion intensities.

Soft soil deposits in Montreal consist mainly of Leda clay, sand, silt and gravel and range in thickness from 1 to 33 m. The shear wave velocity of bedrock is highly variable (from 1000 m/s to 4000 m/s), whereas the V_s profiles of the soft soil deposits vary between 150 m/s and 450 m/s [5]. Consequently, the high variability in the impedance contrast between surface deposits and bedrock can produce larger amplifications than those specified in the NBCC-2015.

Similar studies have been performed in various parts of the world. Aaquib et al. [9] used the 1D site response analysis program DEEPSOIL to calculate the ground motion amplification factor for shallow sites in moderate to low seismicity regions such as Korea, Central and Eastern North America (CENA), and Australia. They noted that peak shear strains occurred at a depth immediately above the interface between the alluvial and weathered soil and that the high impedance contrast at the interface between the alluvial and weathered soil caused the largest amplifications of ground motion. The ratio of maximum amplification factors was shown to generally increase with an increase in the site period [9].

Ciancimino et al. [10] generated a sample of V_s profiles by Monte Carlo simulation and calculated site amplifications with the 1-D Equivalent Linear (EQL) seismic response analysis. The standard deviation of soil amplification factors are estimated as a function of V_s30 and site fundamental frequency for ground motions scaled to PGA 0.1 g. The long period amplification factors do not show significant dispersion (coefficient of variation around 10–20%) and increases with decreasing V_s30 and then decreases when nonlinearity effects become predominant (V_s30 around 200 m/s).

Hashhash et al. [11] developed response spectrum site amplifications by using ten representative V_s profiles for CENA site conditions and a randomization scheme for the shear wave velocity profiles, modulus reduction and damping curves, and rock outcrop ground motions. A high reference rock shear velocity (3000 m/s) was used throughout to account for typical high impedance contrasts in CENA. Both synthetic and rock ground motion records with magnitudes ranging from 4.5 to 7.5 and epicentral distances from 5 to 250 km were used in the simulations.

Harmon et al. [12] used a similar 1-D ground response analysis procedure as in [11] to estimate site amplification functions for CENA and the Western US (WUS). They noted the large stiffness of the bedrock across the CENA. They also noted that the variability in soft soil stratigraphy for a given V_s30 is generally larger for CENA than for WUS sites due to a higher velocity contrast between soft soil layers and bedrock. It is noteworthy that in both papers [11,12] a constant bedrock V_s was used in all the analyses.

Site amplifications in Pistrino di Sotto in Italy were investigated by Sextos et al. [13] to determine how site conditions, such as ancient riverbeds, are related to higher levels of damage observed during a recent earthquake. Rodriquez-Marek et al. [14] used a logic tree approach to evaluate the uncertainties in amplifications. The nodes of the logic trees include shallow V_s profiles, deeper V_s profiles, depth of impedance contrasts, uncertainty on the depth of the rock surface, and uncertainty on the thickness of the weathered rock zone. Falcone et al. [15] analyzed the amplification factors for peak ground acceleration and for the velocity spectrum for a site in Northern Italy (Dovadola). Falcone et al. [16] also performed numerical simulations of seismic site response to produce maps of seismic shaking amplification factors.

Boudghene et al. [17] used an artificial neural network to investigate the correlation between the average amplification factor and various site parameters such as the depth to bedrock, average shear velocity (V_sm) to bedrock, V_s30, shear wave velocity of bedrock (V_bedrock), velocity contrast, and soil fundamental frequency (f₀). They concluded that site fundamental frequency f₀, and V_s30 are the best predictors of site amplification factors F_a and F_v. Zhu et al. [18] also evaluated site effects to obtain the best site parameters in a complement to the V_s30 parameter.

Gobbi et al. [19] performed 1D simulation of site response analysis with a set of 40 earthquake records for stochastically generated horizontally layered soil profiles. They concluded that amplification factors are best predicted as a function of V_s30 in conjunction with the dominant period of the site and the shear wave velocity gradient. Fayjaloun et al. [20] investigated the damage scenario of residential buildings shaken by moderate ground acceleration (200 cm²/s) with and without consideration for the site amplification factor of 1.5. It is noted from Régnier et al. [21,22] that the EQL analysis should be considered reliable with reference to PGA lower than 0.2–0.3 g.

In our study, 1-D EL ground response analyses are performed to investigate the influence of ground motions, nonlinear dynamic soil properties, V_s profiles and impedance contrast on the amplification factors for typical sites in the Greater Montreal area. The results are presented as a function of PGA at rock site for periods of 0.01 s, 0.2 s, 0.5 s and 1.0 s for site classes C and D profiles. The results obtained are compared to the corresponding site factors from NBCC (2015) and from other CENA sites.

2. Methodology of Equivalent Linear 1-D Analysis

Equivalent linear one-dimensional modeling of site response approximates the nonlinear cyclic response of soil through the use of a degradation curve for the shear modulus and a damping ratio curve (e.g., [23,24,25,26,27,28,29,30]). Seismic ground motion propagating through soft soil deposits generates shear waves and Raleigh waves in the subsurface zones of the earth. In practice, shear wave propagation from the bedrock to the ground surface is approximated as 1-D vertically propagating waves from the underlying rock formation.

For layered soil deposits, the first step is to set initial estimates for the shear modulus G and damping for each layer. Then, the transfer function between bedrock and each layer for all major frequencies of the rock motion are computed. Using the product of the Fourier spectra of the bedrock accelerograms and transfer function at various frequencies of the motions, an updated Fourier response spectra is computed for each layer. Next, acceleration, displacement and shear strain time histories are computed for each layer using the inverse fast Fourier technique [31]. After the shear strain time histories are obtained for a layer, an estimate of an effective shear strain (65% of peak strain as proposed in [31]) is determined for the layer. The modulus reduction and damping curves of each layer were then used to obtain the updated values of shear modulus and damping ratio compatible with the current effective strain. Subsequent iterations continue with the updated values until convergence is achieved for the strain levels between two consecutive runs [31].

In our study, it is assumed that the soil layer at each site is underlain by elastic rock. In Montreal, predicted site fundamental frequencies using 1-D EL are in accordance with field measurements [4]. The 1-D Equivalent linear analysis in the frequency domain is used to carry out parametric studies on the effect of ground motion characteristics, dynamic soil properties and impedance contrast ratio for typical soil profiles.

2.1. Selection of Input Rock Ground Motions

A set of 15 records of rock motions, listed in

Table 1. Characteristics of ground input motions selected for the equivalent 1-D analysis in Montreal.

Input Motion	Year	Earthquake Mechanism	Site Class	PGA (g)	Moment Magnitude	Epicentral Distance (km)
Irpinia, Italy	1980	Normal Faulting	B	0.080	6.9	28
Morgan Hil, Gilroy	1984	Strike-slip	A	0.100	6.19	15
Saguenay (Chicoutimi)	1988	Strike-slip	A	0.130	5.7	30
Saguenay (St-Andrea)	1988	Strike-slip	A	0.160	5.7	70
Saguenay (La Malbaie)	1988	Strike-slip	A	0.170	5.7	95
Southern Iran	1990	Strike-slip	B	0.160	6.2	18
Greece	1990	Normal Faulting	B	0.120	5.9	19
Western Iran	1990	Strike-slip	B	0.130	7.4	85
Basilicata, Italy	1998	Normal Faulting	B	0.160	5.6	18
Kocaeli, Turkey	1999	Strike-slip	B	0.26	7.51	11
Tottori, Japan	2000	Strike-slip	B	0.130	6.61	15
Presbytere, Martinique	2007	Normal Faulting	B	0.160	7.4	69
Guadeloupe	2007	Normal Faulting	B	0.080	7.4	144
Goudeloupe, France	2007	Strike-slip	B	0.110	7.4	167
Lipari, Italy	2010	Strike-slip	B	0.32	4.7	12
Synthetic-44	N/A	Strike-slip	A	0.133	7	99
Synthetic-27 *	N/A	Strike-slip	A	0.120	7	70
Synthetic-9 *	N/A	Strike-slip	A	0.280	7	45
Synthetic-14 *	N/A	Strike-slip	A	0.210	7	51
Synthetic-19 *	N/A	Strike-slip	A	0.330	6	27

* the number corresponds to the identification value given to the synthetic accelerogram in [33].

, is selected from the databases of the Pacific Earthquake Engineering Research Center (PEER), the Engineering Strong-Motion database and the Center for Engineering Strong Motion [32]. An additional five synthetic ground motions are selected from [33] that provide a good fit to short periods of the Uniform Hazard Spectra (UHS) for Eastern Canada, which are important for amplification at sites with shallow deposits. The ground motions are selected to match the 2% in 50 years spectrum for Montreal for moment magnitudes between 4.5 and 7.5 and epicentral distances from 10 to 170 km on strike-slip and normal faults recorded at class B and A sites (Figure 1 and Figure 2). The accelerograms are selected to ensure an approximate match to the UHS PGA of 0.25 g, and an overall match with the shape of the UHS for periods from 0.1 s to 2 s.

The recorded ground motions that are selected in Table 1 have a high energy content at short periods, as shown in the spectrum of Figure 1 and Figure 2, which is characteristic for Eastern North America and important for the shallow sites with fundamental frequency between 3 and 20 Hz which are prevalent in Montreal.

The site factors are first computed relative to ground motions on rock obtained at class A and B sites. Since the NBCC-2015 site factors are provided for sites with V_s30 of 760 m/s corresponding to the boundary between soil classes C and B, the site factors for Class D sites calculated in this study are normalized relative class C site factors to be consistent with the NBCC.

2.2. Randomization of V_s Profiles Obtained from Seismic Surveys

The surface geology of Montreal is an interlayered deposit with clay, silt and dense sand overlying till or rock and is the result of several alternating periods of glaciation followed by the emergence of the Champlain Sea and channeling by the St. Lawrence River and its tributaries [34]. Equivalent linear 1-D analyses of ground response are performed for a set of 12 sites representative of soil classes C and D (Figure 3). Their location is shown on Figure 3 with the probabilistic microzonation maps produced by Talukder and Chouinard [35].

The selected sites are representative of the variability in the stratigraphy and depth to bedrock for each site class. Figure 4a,b show the V_s profiles obtained from seismic surveys for 7 sites selected in class C and 5 sites in class D, respectively. The variability of V_s profiles is larger at shallow depths (0 to 10 m) and can be attributed to the heterogeneity of soil deposits compared to deeper layers that consist mainly of clay. The V_s profiles for Montreal are significantly different from the USGS Class C and D profiles [36] and other CENA sites [37] such that sites with similar V_s30 values will have different fundamental frequencies and amplification factors.

Bedrock in Montreal consists mainly of Ordovician limestone, dolomite, shale and Cambrian sandstone and can vary widely in quality [34]. Boyer [38] notes that average V_s for shale and limestone are 967 and 2900 m/s, respectively with high coefficient of variation. The surface geology of Montreal is an interlayered deposit with clay, silt and dense sand overlying till or rock. The unit weight of the soils is listed in

Table 2. Characteristics of soil properties of Montreal sites.

Soil Type	Soil Unit Weight (kN/m3)	Reference
Leda clay	16.5	[30]
Sand	20	[34]
Silt	21	[34]
Till	24	[34]
Bedrock	27.3	[5]

. It is noted from Table 2 that the density of Leda clay is lowest among other soil types.

One of the important sources of variability for amplification factors for a given site class is the variability in soil profiles [39,40,41]. The procedure proposed by Silva et al. [36] is used to generate random soil profiles that account for trends and variability as a function of depth as well as the spatial correlation between shear wave velocities. Shear wave velocities are assumed to be lognormally distributed and interlayer correlation coefficients as a function of depth and layer thickness are proposed from the analysis of over 900 velocity profiles. The soil column is divided in layers and the shear wave velocity for the layer is normalized such that,

$${\mathrm{Z}}_{\mathrm{i}}=\frac{\ln \left({\mathrm{V}}_{\mathrm{i}}\right)-\ln \left({\mathrm{V}}_{\mathrm{median}}\left({\mathrm{h}}_{\mathrm{i}}\right)\right)}{{\mathsf{\sigma }}_{\mathrm{lnV}}}$$

(1)

where, V_i is the shear-wave velocity of the ith layer, V_median(h_i) is the median shear-wave velocity at mid-depth of the layer and σ_lnV is the standard deviation of ln(V_i). The normalized values Z_i are correlated with the layer above (i − 1) using a first order auto-regressive relation:

$${\mathrm{Z}}_{\mathrm{i}}=\mathsf{\rho }·{\mathrm{Z}}_{\mathrm{i}-1}{ \pm \mathsf{\epsilon }}_{\mathrm{i}}·\sqrt{\left(1-{\mathsf{\rho }}^2\right)}; \mathrm{i} > 1$$

(2)

where ε_i is a normal random variable with zero mean and unit standard deviation; and ρ is the interlayer correlation coefficient which is a function of the depth (h) and thickness (t) of the layer.

$$\mathsf{\rho }\left(\mathrm{h}, \mathrm{t}\right)=\left[1-\mathsf{\rho }\left(\mathrm{h}\right)\right]{\mathsf{\rho }}_{\mathrm{t}}\left(\mathrm{t}\right)+\mathsf{\rho }\left(\mathrm{h}\right)$$

(3)

where ρ(h) is a depth dependent interlayer correlation coefficient, and ρ_t(t) is a thickness dependent correlation coefficient:

$$\mathsf{\rho }\left(\mathrm{h}\right){=\mathsf{\rho }}_{200}{\left(\frac{\mathrm{h}}{200}\right)}^{\mathrm{b}}; \mathrm{h} \le 200 \mathrm{m}$$

(4)

where, ρ₂₀₀ is the correlation coefficient at 200 m depth.

The variability (σ_lnVs) of shear wave velocities for each soil class for sites in California is between 0.27 and 0.37. The correlation functions are considered to be valid for other locations and have been used by several other authors for randomizing soil profiles for the Eastern and Central United States [40,41,42]. Since the V_s profiles for Montreal do not exceed 30 m, Equation (4) can be used to calculate the correlation coefficient of the V_s randomization. As a result, the information required to develop randomized shear-wave velocity profiles consists of the baseline shear-wave velocity profiles, the standard deviation of the natural log of the shear-wave velocity σ_lnV, and the interlayer correlation model parameter ρ(h) [36]. The software Strata [8] is used to randomly generate soil profiles using this method for site classes C (Figure 5) and D (Figure 6). A set of 30 V_s profiles were generated for each site class. A total of 210 V_s profiles were generated for the 7 sites in class C (Figure 5) and 150 V_s profiles were generated for the 5 sites in class D (Figure 6). In Figure 5 and Figure 6, the site fundamental frequency T₀ is referred to the target V_s profile. The T₀ for each site is obtained directly from the Microtremor H/V measurements.

2.3. Randomization of G/G_max and Damping Ratios Curves

A randomization procedure is used to account for the uncertainty associated with the dynamic properties of each soil type. In the case of the normalized modulus, the uncertainty is maximum in the middle range since the curve is asymptotic to 0 and 1 at the ends. A normal distribution and the parameters proposed by Darendeli [29] for G/G_max and the damping ratio D are used in this application:

$${\mathsf{\sigma }}_{\mathrm{G}/\mathrm{Gmax}}=0.015+0.16·\sqrt{\left(0.25-{\left(\mathrm{G}/{\mathrm{G}}_{\max }-0.5\right)}^2\right)}$$

(5)

$${\mathsf{\sigma }}_{\mathrm{D}\%}=0.0067+0.78·\sqrt{\mathrm{D}(\%)}$$

(6)

The non-linear properties G/G_max and damping ratios as a function of shear strain γ (%) are negatively correlated and pairs of randomized curves are obtained with the following expressions as a function of shear strain [1]:

$$\left[\begin{array}{c}{\mathrm{G}/\mathrm{G}}_{\max }\left(\mathsf{\gamma }\right)\\ \mathrm{D}\left(\mathsf{\gamma }\right)\end{array}\right]=\left[\begin{array}{cc}{\mathsf{\sigma }}_{\mathrm{G}/\mathrm{Gmax}}& 0\\ {\mathsf{\rho }}_{\mathrm{G}/\mathrm{Gmax},\mathrm{D}}& {\mathsf{\sigma }}_{\mathrm{G}/\mathrm{Gmax}}\sqrt{1-{\left({\mathsf{\rho }}_{\mathrm{G}/\mathrm{Gmax},\mathrm{D}}\right)}^2}\end{array}\right]\left[\begin{array}{c}{\mathsf{\epsilon }}_1\\ {\mathsf{\epsilon }}_2\end{array}\right]+\left[\begin{array}{c}{\mathrm{m}}_{{\mathrm{G}/\mathrm{G}}_{\max }\left(\mathsf{\gamma }\right)}\\ {\mathrm{m}}_{\mathrm{D}\left(\mathsf{\gamma }\right)}\end{array}\right]$$

(7)

where ε₁ and ε₂ are uncorrelated standard normal random variables with ρ_G/Gmax,D = −0.5. Results for G/G_max and the damping ratios are presented respectively in Figure 7 and Figure 8 for clay, sand, silt and gravel.

The shear modulus G/G_max and damping ratio curves for clay [30], sand [29], silt [26] and gravel [24] used in this study are shown in Figure 7 and Figure 8, respectively. Each V_s profile is associated with specific G/Gmax and damping ratio curves that are related to the stiffness of the soil layers. For example, the stiffest V_s profile is associated with the highest G/Gmax curve and the lowest damping ratio curve.

2.4. Relationship between Bedrock V_s versus Rock Quality Designation (RQD)

The rock quality designation is a standard parameter in drill core logging and is a leading indicator for low-quality rock zones [43]. Data on the rock shear wave velocity at various sites where the RQD was also available show a strong correlation between the two measures (Figure 9):

$${\mathrm{Vs}}_{\mathrm{rock}}=816.12+422\ast \mathrm{F}; \mathsf{\sigma }=322 \mathrm{m}/\mathrm{s}; {\mathrm{R}}^2=0.72$$

(8)

where F is the category from 1 to 5 corresponding to the standard five ranges of RQD (0–25%, 25–50%, 50–75%, 75–90%, 90–100%). Rosset et al. [5] considered the average and standard deviation for the rock shear velocity in Montreal of 2300 m/s and 590 m/s, respectively without the RQD information. Here, V_s for rock depends on RQD and the standard deviation is reduced to 322 m/s (Figure 9). This information is used to obtain more precise estimates of site amplifications given the strong influence of the impedance ratio [9,11,12]. Consequently, site amplification factors are obtained in the following sections for RQD classes 1, 3 and 5 (0–25%, 50–75% and 90–100%, respectively). Falcone et al. [44] estimates the effect of stiffness and thickness of soft bedrock (V_s ≤ 800 m/s) on surface ground motion modifications.

3. Site Amplification Sensitivity Analysis

The 1-D ground response of the 12 soil profiles in class C and D with the randomized V_s profiles (Figure 5 and Figure 6) are evaluated by the EL 1D method using Strata [8] under total stress conditions. The set of 20 rock motions are scaled from 0.1 to 0.5 g by increments of 0.1 g. Spectral acceleration response calculations are performed for 5% damping to be consistent with NBCC (2015). The uncertainty associated with the RQD classification is estimated by applying the Rosenblueth point estimation procedure [45,46] in combination with Monte Carlo simulation to reduce the total number of simulations performed. A total of 63,000 and 45,000 simulations have been performed for site classes C and D, respectively.

Figure 10, Figure 11 and Figure 12 show the amplification factor at low period (T = 0.01 s) as a function of V_s30 and site fundamental frequency F₀ for different RQDs of 1, 3 and 5, and PGA of 0.1 g. Amplifications for other PGA values are provided in [47]. They are compared with those obtained previously by Chouinard and Rosset [48] and show good agreement. It is observed from Figure 10, Figure 11 and Figure 12 that for all 3 RQDs, a peak is observed in the non-linear relation between F(0.01 s) and V_s30 when V_s30 is in the range of 360–550 m/s.

Tempa et al. [49] conducted a parametric study with 1-D EL method to investigate how the bedrock characteristics influence the surface ground motions for various depths of soil profiles in Bhutan. The PGA values near ground surface are in the range between 0.191 and 0.221 g for 30 m thick soil profiles [49]. To identify the amplification characteristics due to geotechnical uncertainty in site classification for South Korean sites, Sun et al. [50] proposed a dual framework: firstly, a site classification method by using depth to bedrock, V_s30 and site period, secondly, a site classification method by using digital elevation, SPT-N values and slope of surface soil (topography). Kim et al. [51] investigated the correlation between the earthquake damage mechanism and site effects by evaluating spatial comparison earthquake damage and seismic hazard. GIS-based method consisting of borehole logs, a digital elevation model (DEM), an administrative division map, and geological maps was used by Kim et al. [52] in order to present site-specific surface and subsurface conditions in Daejeon city. Chieffo et al. [53] evaluated the influence of site amplification on the seismic vulnerability of historical structure in the Molise Region of Italy by analyzing distributions of earthquake intensities and seismogenic-faults.

Confidence intervals are provided for the mean amplification factor F(0.01 s) as well as for predicted amplification factor. The latter corresponds to the uncertainty expected for single events and highlights the possibility of amplifications significantly larger than mean amplifications specified in design codes for specific earthquakes.

Figure 13, Figure 14 and Figure 15 show the amplification factors at the periods of 0.2 s, 0.5 s and 1.0 s as a function of V_s30 for RQDs of 1, 3, and 5, and PGA of 0.1 g. Amplifications for PGAs from 0.2 g to 0.5 g are provided in Talukder [47]. The amplification factors for the structural periods of 0.2 s, 0.5 s and 1.0 s increase as a function of RQD, which is consistent with an increase in the impedance contrast [9,11,12,54]. Site amplification factors are compared with those recommended in NBCC (2015) values and those obtained previously by Kim and Yoon [54]. In the NBCC (2015), site factors are normalized relative to site class C amplifications (F_C = 1).

For comparison purposes, the NBCC site factors for class A site (F_A) and class B site (F_B) have been normalized relative to rock amplifications (F_A for RQD of 3 and 5 and F_B for RQD of 1). It is noted from Figure 13a, Figure 14a and Figure 15a that the results show higher amplifications than those of [54] and NBCC 2015 for a short period (T = 0.2 s), especially for high RQD categories. The results show amplifications similar to those of [54] and NBCC (2015) for longer periods (T = 0.5 s and 1.0 s) for all RQD categories (Figure 13a,b, Figure 14a,b and Figure 15a,b).

Figure 16 compares median amplifications at periods of 0.01 s, 0.2 s, 0.5 s and 1.0 s for Montreal to CENA median site factors [37]. The CENA site factors are obtained for soil deposits of 15 and 33 m, site classes C and D, and bedrock velocity of 2900 m/s [37], which corresponds to a RQD of 5 in Montreal. The results show good agreement with the results for site class C for low periods of 0.01 s and 0.2 s (15 m deposits are typical for site class C sites in Montreal). For class D sites, the amplifications for Montreal are higher than those for CENA and the agreement is closer for high periods and 33 m deposits, which are typical for class D sites in Montreal. Higher amplifications can be explained by the low damping ratio of Leda clay in comparison to other soft soils.

Non-linear effects are more important for sites with low V_s30 and smaller amplifications can be observed, specially for larger PGA values (Figure 16). It can be noted from Figure 16c,d that for softer soils (class D sites), higher amplification factors are obtained at longer periods of T = 0.5 s and 1.0 s. By contrast, for stiffer soils (class C sites), relatively lower amplification factors are obtained for shorter periods of T = 0.01 s and 0.2 s. Similar variation of the amplification factors as a function soil site class are reported in recent literature [44,55,56,57]. Akkar et al. [55] formulated new ground motion prediction equations for Europe and Middle East by including non-linear site amplification function that is a function of V_s30 values from 150 to 1200 m/s and reference peak ground acceleration. Bindi et al. [56] observed it from their 1-D ground response analyses that site class A (rock) dominates the site-to-site variability for short site fundamental periods while for long periods the difference among the classes is less pronounced, with soil sites (classes C and D). Kotha et al. [57] presented nonlinear curves site amplification factors for periods of T = 0.01, 0.1, 1 s as a function of with V_s30. They observed that short-period site response at rock sites shows high variability compared to softer soils but the long-period site-response of rock sites show less variability than that of softer soils (e.g., 180 m/s ≤ V_s30 < 800 m/s).

The amplification factors F_D and F_C corresponding to V_s30 = 250 m/s and 450 m/s are used to define representative values for soil classes D and C, respectively, in NBCC (2015). The amplification factors F_D and F_C are assumed to be log-normally distributed and uncorrelated.

$$\left(\begin{array}{c}{\mathrm{lnF}}_{\mathrm{D}}\\ {\mathrm{lnF}}_{\mathrm{C}}\end{array}\right) ~ \mathrm{N}\left[\begin{array}{cc}\left(\begin{array}{c}{\mathsf{\mu }}_{{\mathrm{lnF}}_{\mathrm{D}}}\\ {\mathsf{\mu }}_{{\mathrm{lnF}}_{\mathrm{C}}}\end{array}\right),& \left(\begin{array}{cc}{\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{D}}}^2& 0\\ 0& {\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{C}}}^2\end{array}\right)\end{array}\right]$$

(9)

where,

$${\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{C}}}^2=\ln \left[1+\frac{{\mathsf{\sigma }}_{{\mathrm{F}}_{\mathrm{C}}}^2}{{\mathsf{\mu }}_{{\mathrm{F}}_{\mathrm{C}}}^2}\right]; {\mathsf{\mu }}_{{\mathrm{lnF}}_{\mathrm{C}}}=\ln \left[{\mathsf{\mu }}_{{\mathrm{F}}_{\mathrm{C}}}\right]-0.5{\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{C}}}^2$$

$${\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{D}}}^2=\ln \left[1+\frac{{\mathsf{\sigma }}_{{\mathrm{F}}_{\mathrm{D}}}^2}{{\mathsf{\mu }}_{{\mathrm{F}}_{\mathrm{D}}}^2}\right]; {\mathsf{\mu }}_{{\mathrm{lnF}}_{\mathrm{D}}}=\ln \left[{\mathsf{\mu }}_{{\mathrm{F}}_{\mathrm{D}}}\right]-0.5{\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{D}}}^2$$

The mean value of the ratios F_D/F_C, F_A/F_C and F_B/F_C are used in NBCC (2015) to normalize the short period and long period site factors (formerly F_a and F_v), which are also lognormally distributed. The mean value and variance of the ratios can be derived from the distributions of F_D and F_C as:

$${\mathsf{\mu }}_{{\mathrm{F}}_{\mathrm{D}}/{\mathrm{F}}_{\mathrm{C}}}=\exp \left({\mathsf{\mu }}_{{\mathrm{lnF}}_{\mathrm{D}}}-{\mathsf{\mu }}_{{\mathrm{lnF}}_{\mathrm{C}}}+0.5\left({\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{D}}}^2{+\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{C}}}^2\right)\right)$$

(10)

$${\mathsf{\sigma }}_{{\mathrm{F}}_{\mathrm{D}}/{\mathrm{F}}_{\mathrm{C}}}^2=\exp \left(2·\left({\mathsf{\mu }}_{{\mathrm{lnF}}_{\mathrm{D}}}-{\mathsf{\mu }}_{{\mathrm{lnF}}_{\mathrm{C}}}\right){+\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{D}}}^2{+\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{C}}}^2\right)·\left(\exp \left({\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{D}}}^2{+\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{C}}}^2\right)-1\right)$$

(11)

and the 95% confidence interval on F_D/F_C as:

$$\exp \left[\left({\mathsf{\mu }}_{{\mathrm{lnF}}_{\mathrm{D}}}-{\mathsf{\mu }}_{{\mathrm{lnF}}_{\mathrm{C}}}\right) \pm 1.96\sqrt{\left({\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{D}}}^2{+\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{C}}}^2\right)}\right]$$

(12)

A feature of amplifications obtained by simulation as opposed to observed data is that each scaled ground motion record is applied to a set of randomized soil profiles and dynamic soil properties. This introduces positive correlation between amplifications obtained for different soil classes that are not present in field observations of amplifications. This correlation is a result of using the same input ground motion for each randomized profile. Considering the correlations, the mean value and the variance for the amplifications normalized relative to site class are as follows:

$$\left(\begin{array}{c}{\mathrm{lnF}}_{\mathrm{D}}\\ {\mathrm{lnF}}_{\mathrm{C}}\end{array}\right)~\mathrm{N}\left[\begin{array}{cc}\left(\begin{array}{c}{\mathsf{\mu }}_{{\mathrm{lnF}}_{\mathrm{D}}}\\ {\mathsf{\mu }}_{{\mathrm{lnF}}_{\mathrm{C}}}\end{array}\right),& \left(\begin{array}{cc}{\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{D}}}^2& {\mathsf{\rho }\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{D}}}{\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{C}}}\\ {\mathsf{\rho }\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{D}}}{\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{C}}}& {\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{C}}}^2\end{array}\right)\end{array}\right]$$

(13)

$${\mathsf{\mu }}_{{\mathrm{F}}_{\mathrm{D}}/{\mathrm{F}}_{\mathrm{C}}}=\exp \left({\mathsf{\mu }}_{{\mathrm{lnF}}_{\mathrm{D}} }-{ \mathsf{\mu }}_{{\mathrm{lnF}}_{\mathrm{C}}}+0.5\left({\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{D}}}^2-{2\mathsf{\rho }\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{D}}}{\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{C}}}{+\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{C}}}^2\right)\right)$$

(14)

$${\mathsf{\sigma }}_{{\mathrm{F}}_{\mathrm{D}}/{\mathrm{F}}_{\mathrm{C}} }^2=\exp \left(2·\left({\mathsf{\mu }}_{{\mathrm{lnF}}_{\mathrm{D}}}-{ \mathsf{\mu }}_{{\mathrm{lnF}}_{\mathrm{C}}}\right){+\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{D}}}^2{+\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{C}}}^2-{2\mathsf{\rho }\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{D}}}{\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{C}}}\right)·\left(\exp \left({\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{D}}}^2{+\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{C}}}^2-{2\mathsf{\rho }\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{D}}}{\mathsf{\sigma }}_{{\mathrm{lnF}}_{\mathrm{C}}}\right)-1\right)$$

(15)

where ρ is the correlation coefficient between lnF_D and lnF_C. The coefficient of correlation coefficient as a function of PGA for RQDs of 1, 3 and 5 and a period of T = 0.01 s, 0.2 s, 0.5 s, 1.0 s are shown in Figure 17.

The correlation coefficient is not very significant for a period of 0.01 s (Figure 17a) corresponding to RQD = 1 and decreases as a function of PGA. The correlation coefficient is more significant for longer periods of 0.5 s and 1.0 s (Figure 17b) and tends to decrease with increasing PGA. The correlation is largest for higher periods due to similarities in soil responses for the two site classes at small PGA, decreases with increasing ground motions due to non-linear effects. Figure 18 presents the site factor E(F_D/F_C) for a period of 0.01 s as a function of PGA as well as the 95% confidence interval. Figure 19, Figure 20 and Figure 21 compare the site factor for RQD of 1, 3 and 5 with the 2015 NBCC site factors for periods of 0.2 s, 0.5 s, and 1.0 s.

The ratio of E(F_D/F_C) for T = 0.01 s is slightly less than unity and is generally in agreement with the 2015 NBCC site factors (Figure 18). The amplification factors at periods of 0.01 s (Figure 18) and 0.2 s (Figure 19) are of the same order of magnitude as those of the NBCC 2015 with a slight positive trend as a function of PGA. The differences are much more significant for periods of 0.5 s (Figure 20) and 1.0 s (Figure 21) and indicate that NBCC 2015 underestimates amplifications in this range. The 95% confidence intervals on the amplification factors are fairly wide and indicate that there is a high probability of much larger amplifications that those specified in design codes, which could increase the level of damage and expected losses in seismic risk analyses given the asymmetry of damage as a function of amplification.

4. Discussion

Average site amplification factors for PGA of 0.1 g are presented in order to characterize linear site amplifications. For this level of ground motions, the results indicate that the average amplifications increase with RQD (or the impedance contrast). Another feature of the results for a period of 0.2 s (Figure 13a, Figure 14a and Figure 15a) is that the average amplifications do not monotonically decrease as a function of V_s30 and that for some class D sites; the average amplifications can be smaller than for class C sites. Similar results have been reported in the recent literature of NGA-East [58,59].

Figure 13, Figure 14 and Figure 15 compare the calculated site amplification factors as a function of V_s30 to the site factors of NBCC (2015) and with the mean amplifications from [54] for Korean sites. Results indicate again that NBCC (2015) underestimates the mean amplifications and that the best agreement corresponds to a RQD of 1 or lower impedance ratios. This result is consistent with the general observation that impedance ratios are smaller for western sites that form the basis for NBCC (2015). Good agreement is obtained with the results of Kim and Yoon [54]. The calculated site factors for T = 0.2 s, 0.5 s and 1.0 s for ground motion records scaled to 0.2, 0.3, 0.4 and 0.5 g show similar tendencies but with lower site factors and better agreement with NBCC due to non-linear effects.

As presented in Figure 16, the comparison of the results for RQD = 5 with those of [37] for shallow soil class C and D sites (15 m and 33 m) in CENA shows comparable trends for short period (T = 0.01 s and 0.2 s). Amplification factors for class C sites are slightly higher than for class D sites and the effect of the deposit thickness has a limited influence for the range of depths considered. For longer periods (T = 0.5 s and 1.0 s), amplification factors for class D sites are higher than those for class C sites and an increase in the thickness of the deposits increases the amplification. Results are shown for a RQD of 5 since the rock velocities for this category is 2900 m/s, which is similar to the velocity assumed in the CENA study (3000 m/s). CENA results are slightly smaller than those obtained with the current study, which could be partially explained by the lower degree of damping for Leda clay and differences in shear wave velocity profiles in the CENA site profiles (Figure 4). The soil profile for class D sites in the CENA study (Figure 4a) is similar to the median profile for sites in Montreal. Conversely, the soil profile for class C sites in the CENA study (Figure 4b) is stiffer than all of the class C site profiles used in this study, which explains partially the higher amplification factors obtained due to a higher impedance ratio.

The comparison of normalized amplification factors of F_D/F_C (Figure 18, Figure 19, Figure 20 and Figure 21) with NBCC (2015) indicates that the normalized factors are similar for short periods (T = 0.01 s and 0.2 s) and significantly larger for longer periods (T = 0.5 s and 1.0 s). The uncertainty on the ratio of F_D/F_C can be large due to uncertainties in ground motions, soil profiles and non-linear dynamic properties, which induce a wide range of site amplification values for individual earthquake occurrences.

The normalized amplification factors of F_D/F_C (Figure 19, Figure 20 and Figure 21) are also compared with the site factors of Hashash et al. [11]. For PGA from 0.1 to 0.5 g and for RQD of 1 to 5, the normalized factors calculated at T = 0.2 s (Figure 19) are greater than the site factors calculated by [11] at T = 0.1 s but smaller for the short period T = 0.3 s. The normalized mean factor F_D/F_C calculated at T = 0.2 s for PGA of 0.1 g is significantly smaller than the one in [11] for T = 0.3 s. For the range of rock PGA of 0.15 to 0.5 g, the site factors in [11] are found to be within the 95% prediction interval calculated for the mean F_D/F_C. The mean F_D/F_C at T = 0.5 s (Figure 20) and the ones in [11] for periods of T = 0.3 s and 1.0 s are similar with increasing PGA, except for PGA values less than 0.2 g. It is observed from Figure 21 that the mean factor F_D/F_C at T = 1.0 s is similar to the site factor in [11] for the same period. Figure 20 and Figure 21 show that the site factors calculated in [11] are within the 95% prediction interval of the mean F_D/F_C for PGA between 0.1 and 0.5 g. Major sources of difference between the amplifications calculated in this study and those by Hashash et al. [11] are the variation of impedance contrast between soil and rock interface, the frequency content of input rock motions, and V_s profiles.

5. Conclusions

The 1-D EL modeling of the propagation of seismic ground motions was performed for a collection of soil profiles that are classified as site classes C and D using the V_s30 criterion. Uncertainty in site responses was investigated by randomizing each of the soil profiles, bedrock properties, non-linear dynamic properties and by selecting a set of representative ground motions. Results were summarized such that they can be compared to currently site factors specified in the NBCC (2015) and other published results. For this purpose, amplifications were obtained for periods of 0.01 s, 0.2 s, 0.5 s and 1.0 s and for ground motions scaled at 0.1 to 0.5 g by increment of 0.1 g. The results were presented for site classes based on the V_s30 classification. Data on seismic velocities in bedrock for Montreal indicate a strong correlation between the degree of defects in the rock (measured by the rock quality designation) and shear wave velocity through rock, which affects the impedance contrast. To analyze this effect, results were obtained for the full range of RQD categories (i.e., 1, 3 and 5).

The results indicate that average site amplification factors based on the V_s30 classification are underestimated in the NBCC (2015). This finding is consistent with results obtained for sites along the CENA (e.g., South Carolina) and other regions of the world (e.g., South Korea) that exhibit similar characteristics (shallow soil layers on bedrock of high V_s) of soil layers. The calculated site amplification factors also indicate that there are large uncertainties associated with the site amplifications, which are also consistent with the degree of uncertainty reported for empirically derived amplifications from seismic observations. These uncertainties are typically neglected when performing seismic hazard analyses which may underestimate expected and extreme levels of damage.

Non-linearity in the analysis of Equivalent Linear 1-D site response is taken into account by introducing the nonlinear G/Gmax and damping ratios curves. In this method, it is assumed that the shear strain compatible Shear Modulus and Damping ratio values remains constant throughout the duration of the seismic excitation. This assumption is not fully applicable to a case when loose saturated soil profiles undergo heavy shaking (PGA > 0.3 g). A complete loss of effective stress may occur in loose saturated deposit of sand during earthquakes with PGA > 0.3 g. For earthquakes with PGA >0.3 g, the adoption of a stress-strain hysteretic nonlinear (NL) model in time domain may be considered. It is to be noted that cohesive soils (clay and silt) at Montreal sites are stiff and cohesionless soils (sand and gravel) are very dense (Table 2). In this study, all simulations with input motion PGA >0.3 g have been performed by using the EL method instead of the NL method because the preparation of input data to the EL method requires only experimental G/Gmax and damping ratio curves that are obtained from literature. On the other hand, the parameters that describe the NL models could not be obtained because of the unavailability of substantial field and laboratory data on the Montreal sites. This makes the use of the EL method of analysis for PGA >0.3 g more practical.

The most significant finding of the study is that the published site factors for CENA are shown in this study to be within the 95% prediction interval calculated for the normalized amplification factors of F_D/F_C. It is shown in the study that positive correlation between lnF_D and lnF_C due to ground motion and that the 95% prediction intervals on the normalized mean amplification that are calculated with correlation coefficients for PGA of 0.1 to 0.5 g are seen to be smaller than the 95% prediction intervals calculated without the correlation coefficients of F_D/F_C.