Vs30 Mapping of the Greater Montreal Region Using Multiple Data Sources

The metropolitan community of Montreal (MMC) is located in Eastern Canada and included in the western Quebec seismic zone characterized by shallow crustal earthquakes and moderate seismicity. Most of the urbanized areas are settled close to the Saint-Lawrence River and its tributaries and within the region, delimiting the extension of the clay deposits from the Champlain Sea. The influence of these recent and soft deposits on seismic waves has been observed after the 1988 M5.8 Saguenay earthquake and has proven to be crucial in seismic hazard analysis. The shear-wave velocity Vs averaged over the 30 m of soil, abbreviated Vs30, is one of the most used parameters to characterize the site condition and its influence on seismic waves. Since 2000, a site condition model has been developed for the municipalities of Montreal and Laval, combining seismic and borehole data for risk mitigation purposes. The paper presents an extended version of the Vs30 mapping for the entire region of the MMC, which accounts for half of the population of Quebec, including additional ambient noise recordings, recently updated borehole datasets, geological vector map and unpublished seismic refraction data to derive Vs profiles. The estimated Vs30 values for thousands of sites are then interpolated on a regular grid of 0.01 degrees using the inverse distance weighted interpolation approach. Regions with the lowest estimated Vs30 values where site amplification could be expected on seismic waves are in the Northeastern part and in the Southwest of the MMC. The map expresses in terms of site classes is compared with intensity values derived from citizen observations after recent felt. In general, the highest reported intensity values are found in regions with the lowest Vs30 values on the map. Areas where this rule does not apply, should be investigated further. This site condition model can be used in seismic hazard and risk analysis.


Introduction
The urbanized area of Montreal (Canada) is located in Southeastern Canada, a region with a moderate seismic hazard [1,2]. The city experienced its largest historical earthquake on September 16, 1732 (M5.8), when 300 houses were reported heavily damaged [3]. At that time, the population was only 3000 inhabitants and houses were all located in the old port area. The current population of the Montreal Metropolitan Community (MMC) is around 4 million in an urbanized area of more than 4370 km 2 . In the event of a repetition of the 1732 earthquake, Rosset et al. [4] have estimated that 12% of the residential buildings could expect extensive to complete damage. The repair and replacement costs could reach 12% of the total portfolio value in Montreal and 0.04% in surrounding municipalities, with nonstructural damage accounting for 80% of the total damage on average. The national seismic risk model for Canada recently proposed by Hobbs et al. [5] concludes that Montreal is the first city at risk in Canada after Vancouver. In both analyses, it is noted that the distribution of the soil conditions partially influences the distribution of damage, especially in regions with glacial deposits and recent soft soil deposition, such as the MMC.
Near-surface soft soil deposits, such as clay or sand, can amplify or de-amplify seismic waves, and the velocity contrast between two layers can result in amplification and shortening of the shear-wave wavelength as it passes from one layer to another.
The spatial resolution of the seismic hazard and risk study can be a factor in determining the level of detail required for the site characterization. Studies performed at the national or regional scale have often relied on proxy models, such as the USGS topographic slope-based site classification [6]. However, such models have been shown to perform poorly in regions of past glaciation, such as the Montreal region (e.g., [7,8]). In the latter case, several methods of surface wave analysis have been developed to estimate site-specific shear-wave velocity profiles or V s30 (e.g., guidelines by Foti et al. [9]). The choice of methods is a function of the desired level of detail, time and budget.
In the region of Montreal, the surface geology is composed of sand from the Saint-Lawrence River, marine clay from the Champlain Sea, glacial tills from three separate glaciation periods and bedrock (limestone and shale) from the Cambrian-Ordovician period. The thickness of post-glacial deposits is highly variable within the MMC and can reach depths up to 50 m in the Northeastern part of the region. The Leda clay originates from the Champlain Sea and is very sensitive to ground motions [10]. Amplification effects of the clay have been observed during the M5.8 1988 Saguenay earthquake and from the analysis of felt reports for recent earthquakes near Montreal [11].
In 2001, a seismic microzonation project was initiated in the Montreal region by compiling data from boreholes, seismic surveys, and ambient noise recordings analyzed using the Horizontal-to-Vertical Spectral Ratio (HVSR) method [12,13]. In Canada, a similar approach was followed by Motazedian et al. [14] for the Ottawa region using shear-wave seismic data in 15 cased boreholes, 686 sites investigated using the refraction-reflection seismic method, and recorded ambient noise at about 400 sites when soft-soil thickness exceeded 10 m. Of these, 185 were collocated in sites investigated using the former methods. Over 25 km of continuous seismic sections with V s30 measurements were obtained using the 24-48 channel streamer of horizontal geophones mounted on sleds behind a horizontally polarized swept-frequency vibrator source (Minivib). Recordings from 16 local small earthquakes at soil and bedrock stations complemented the field measurements. A 3D geological model was constructed from the three surficial geological units (clays, sand and till) identified in each of the 21,000 boreholes in the database. Molnar et al. [15] have a similar project in Vancouver based on field measurements and the development of a high-resolution 3D model for later use in numerical simulations of wave propagation. In Europe, Oliveira et al. [16] took a different approach to the microzonation of Lisbon by first compiling a dataset of more than 8000 boreholes with N30 (the average Standard Penetration Test for the first 30 m) as a proxy for V S30 . The 1:10,000 scale geological map was then used to refine the map, which was finally validated by non-invasive field surveys (surface wave seismic and HVSR methods).
These microzonation are based on V S30 , which is the average shear-wave velocity from the ground surface to a depth of 30 m. The calculation of V S30 is as follows: where Z i is the thickness and averaged V si is the shear-wave velocity over the layer i. V Srock is Vs for the rock if it is present in the first 30 m of the profile. The National Building Code of Canada [17] groups this parameter into site classes corresponding to different soil types, as shown in Table 1. The approach, first applied on the island of Montreal (Figure 1), is being extended to the region of the MMC with the benefit of open access data from public institutions in Quebec and newly investigated sites with ambient noise records. Section 2 describes the past and new data collected to perform the mapping in the MMC, as Section 3 details the procedures adopted to calculate V s30 , which depends on the source of collected data. The results of their interpolated values using the inverse distance weighted method are proposed considering a similar weight for the different sources of data. The resulting V s30 maps are discussed in regard to the newly acquired felt reports from citizens after weak earthquakes in the region. The resulting map is used to incorporate the site conditions into the seismic hazard that is calculated from the national probabilistic models or scenario earthquakes.
Geosciences 2023, 13, 256 3 of 16 Quebec and newly investigated sites with ambient noise records. Section 2 describes the past and new data collected to perform the mapping in the MMC, as Section 3 details the procedures adopted to calculate Vs30, which depends on the source of collected data. The results of their interpolated values using the inverse distance weighted method are proposed considering a similar weight for the different sources of data. The resulting Vs30 maps are discussed in regard to the newly acquired felt reports from citizens after weak earthquakes in the region. The resulting map is used to incorporate the site conditions into the seismic hazard that is calculated from the national probabilistic models or scenario earthquakes.

Seismic Measurements
Over the last 20 years, 2653 sites were investigated with ambient noise records in order to estimate the predominant resonance frequency estimated using the Horizontal to Vertical Spectral Ratio (HVSR) method [18]. The locations were typically selected on the

Seismic Measurements
Over the last 20 years, 2653 sites were investigated with ambient noise records in order to estimate the predominant resonance frequency estimated using the Horizontal to Vertical Spectral Ratio (HVSR) method [18]. The locations were typically selected on the basis of the surface geology and favored areas with clay, peat and sand deposits with fewer measurements in zones with tills and rock outcrops. Three different configurations of 24-bit digitizers connected to 3-component velocimeters were used: ORION from Nanometrics Ltd. (Kanata, ONT, Canada). Connected to a Guralp CMG-40T/30 s, CityShark from LEAS connected to a LENNARTZ LE-3D/5s, and Tromino ® stations. A clear peak could be interpreted as the predominant resonance frequency for about twothirds of the sites investigated using the SESAME methodology [19]. The analysis for the remaining sites was more complex and required identifying the frequency peak corresponding to the thickness of soft layers at the closest borehole sites. Rosset et al. [12] described the acquisition procedure and the adopted parameters for the analysis. The sites in the municipalities outside Montreal, Laval and Longueuil were predominantly selected in urbanized zones with a thicker layer of clay based on existing borehole data (as explained in Section 2.2). A set of 13 High-Resolution Multichannel Seismic Reflection profiles using the Minivib System, a mobile "landstreamer" consisting of 48 receivers on small metal sleds (HRSR, see for more details [20]) was done in collaboration with Natural Resources Canada to estimate shear-wave velocity Vs. Multi Acquisition Seismic Waves (MASW) surveys were performed at 29 sites with an acquisition system connected to 24 vertical low-frequency (4.5 Hz) geophones deployed on a line at regular spacing. Vertical stacking of eight impacts with a sledgehammer on a metallic plate was used to generate a seismic image, including Rayleigh waves propagating at the surface. A minimum of two active MASW field records were collected at each site with different source offsets and receiver spacing. Passive recordings were also used to improve the analysis at low frequencies and the dispersion image [21]. Finally, two sites with boreholes were selected to calculate Vs profiles using the downhole seismic method described by Hunter and Crow [22]. The sites were selected because of the availability of cased boreholes at the time of the field campaigns (two of them are distant less than 300 m). The map of Figure 1 locates the investigated sites with different methods.

Borehole and Geological Data
A borehole and geological data have been collected to complement the direct and indirect Vs measurements. The map in Figure 2 shows the location of the data grouped by sources as follows: -Our own dataset was built from borehole reports collected in the cities of Montréal, Laval and Longueuil. - The dataset of the Ministère des Transports Québec (MTQ) includes P-wave seismic refraction (SR) profiles in 1584 sites within the MMC [23]. Each profile provides with thickness and depth of the soil layers. 36% of the sites have a profile with a layer indicated as "unknown type", which overlays a clay layer or another layer with an unknown type of a higher thickness and a different Vp value. The conversion of Vp to S-wave velocity Vs is performed in order to calculate V s30 using Equation (1) (see The MTQ dataset also includes thousands of digitized borehole reports geographically located [23]. - The dataset of the hydrological information system from the Ministère de l'Environnement et de la Lutte contre les Changements Climatiques (MELCC) includes digitized wells reports, as well as a database providing the thickness and soil types along the borehole (donneesquebec.ca). - The location of the outcropping bedrock from the digital version of the 1:2,000,000 geological map of Quebec was updated in 2022 [24].

Vs Derived from the Seismic Refraction Data
The SR data provide Vp profiles which need to be converted to Vs For that, the database from Hunter J. (unpublished, 2011), which compiles Vp and Vs velocities for different soil layers, is used. It includes data measured at different depths in several sites in Ontario for Leda clay, Pleistocene till (2.6 million to 11,700 years ago), as well as sand and silt of Pleistocene and Holocene (11,700 years to present) ages. These data are used to estimate a Vp/Vs relation with depth for each available soil type, as plotted in Figure 3.
The observations are fitted using a relation of the form: where H is the depth in m. Table 2 lists the best-fitted parameters a and b, the standard deviation around the mean value, as well as the regression coefficient R 2 for each soil type. The dataset for Holocene sand is not included in the table because it contains only four measurements, as shown in Figure 3.

Vs Derived from the Seismic Refraction Data
The SR data provide Vp profiles which need to be converted to Vs For that, the database from Hunter J. (unpublished, 2011), which compiles Vp and Vs velocities for different soil layers, is used. It includes data measured at different depths in several sites in Ontario for Leda clay, Pleistocene till (2.6 million to 11,700 years ago), as well as sand and silt of Pleistocene and Holocene (11,700 years to present) ages. These data are used to estimate a V p /V s relation with depth for each available soil type, as plotted in Figure 3.
The observations are fitted using a relation of the form: where H is the depth in m. Table 2 lists the best-fitted parameters a and b, the standard deviation around the mean value, as well as the regression coefficient R 2 for each soil type. The dataset for Holocene sand is not included in the table because it contains only four measurements, as shown in Figure 3. Geosciences 2023, 13, 256 6 of 16

Vs for Rock Using the Rock Quality Designation (RQD)
The Rock Quality Designation (RQD) is a parameter defining the degree of jointing in the rock and ranging from zero to 100%, the latter being the best quality rock. In order to consider this information when estimating Vs for rock, Talukder et al. [25] proposed a relationship as follows: F is a value given for a range of RDQ percentage (1 for 0-20%, 2 for 20-40%, 3 for 40-60%, 4 for 60-80% and 5 for values higher than 80%). The average value of three is assumed when the RQD parameter is missing in borehole data (Vs30 = 2082 m/s).

Vs for Rock Using the Rock Quality Designation (RQD)
The Rock Quality Designation (RQD) is a parameter defining the degree of jointing in the rock and ranging from zero to 100%, the latter being the best quality rock. In order to consider this information when estimating Vs for rock, Talukder et al. [25] proposed a relationship as follows: F is a value given for a range of RDQ percentage (1 for 0-20%, 2 for 20-40%, 3 for 40-60%, 4 for 60-80% and 5 for values higher than 80%). The average value of three is assumed when the RQD parameter is missing in borehole data (V s30 = 2082 m/s).

Estimation of V s30 from the Different Sources of Data
The proposed microzonation classifies local site response in site classes based on V S30 , which is the averaged shear-wave velocity for the first 30 m of soil as calculated in Equation (1). It follows the classification used in the National Building Code of Canada (NBCC, [17]) and considers five classes from A to E defined by the V s30 value (Table 1). Figure 4 is a flowchart describing the procedure used to calculate V s30 for the different sources of data that are further detailed in this section.

Estimation of Vs30 from the Different Sources of Data
The proposed microzonation classifies local site response in site classes based on VS30, which is the averaged shear-wave velocity for the first 30 m of soil as calculated in Equation (1). It follows the classification used in the National Building Code of Canada (NBCC, [17]) and considers five classes from A to E defined by the Vs30 value (Table 1). Figure 4 is a flowchart describing the procedure used to calculate Vs30 for the different sources of data that are further detailed in this section.  Table 2 gives the parameters of the Equation (2).

Vs30 Estimates from Seismic Refraction Data
Vs30 is calculated using Equation (1) after converting Vp profiles to Vs using the relations Vp/Vs as described in Table 2 and Equation (2) for the rock. In several SR profiles, the soil type in one or two layers is undetermined (indicated as an unknown layer). Wherever possible, the unknown soil types were replaced with those identified in boreholes at similar depths and less than 250 m from the investigated site (a distance that corresponds to the maximum spatial correlation of Vs30). If no other sites are within 250 m, backfill material was assumed when Vp is lower than 750 m/s and sand when Vp is between 750 and 2000 m/s. The unknown layer with a Vp greater than 2000 m/s was assumed to be till, as this is the only type with high Vp values. A fixed Vs value of 155 m/s was assumed for backfill, as proposed by Rosset et al. [12].
The large variability of Vp/Vs values for the Leda clay at shallow depths (5-10 m) in Figure 3a is largely influenced by the level of the groundwater table at the site. As this information is not available and is highly variable in time and space, the mean Vp/Vs value was chosen for all sites, as in Table 2. The sites with a clay layer in the upper 10 m represent 10% of the number of sites having a first layer as clay (which is less than 2% of the total  Table 2 gives the parameters of the Equation (2).

V s30 Estimates from Seismic Refraction Data
V s30 is calculated using Equation (1) after converting Vp profiles to Vs using the relations Vp/Vs as described in Table 2 and Equation (2) for the rock. In several SR profiles, the soil type in one or two layers is undetermined (indicated as an unknown layer). Wherever possible, the unknown soil types were replaced with those identified in boreholes at similar depths and less than 250 m from the investigated site (a distance that corresponds to the maximum spatial correlation of V s30 ). If no other sites are within 250 m, backfill material was assumed when Vp is lower than 750 m/s and sand when Vp is between 750 and 2000 m/s. The unknown layer with a Vp greater than 2000 m/s was assumed to be till, as this is the only type with high Vp values. A fixed Vs value of 155 m/s was assumed for backfill, as proposed by Rosset et al. [12].
The large variability of Vp/Vs values for the Leda clay at shallow depths (5-10 m) in Figure 3a is largely influenced by the level of the groundwater table at the site. As this information is not available and is highly variable in time and space, the mean Vp/Vs value was chosen for all sites, as in Table 2. The sites with a clay layer in the upper 10 m represent 10% of the number of sites having a first layer as clay (which is less than 2% of the total number of SR sites). Vs varies between 60 and 100 m/s, values measured for clay at this depth by Mayne et al. [10] at the Canadian test site in Gloucester, Ontario. In total, V S30 is calculated for 1014 SR sites and classified into soil classes based on their values, as shown in the map of Figure 5.

Vs30 Estimates from Borehole Data
The two datasets of the MTQ and the MELCC are used to estimate VS30. Each selected borehole site includes information on the soil types, their relative thicknesses, as well as the RQD of bedrock. Vs values of the four main soil types (backfill, clay, sand and till) were estimated using the proposed Vs relation for the four prevalent soil layers in Montreal provided by Rosset et al. [12]. For clay and sand types, Vs profile with depth H is proposed combining various seismic measurements and borehole data. Table 3 lists the model coefficients, and Figure 6 show the data and model for sand and clay up to 35 m. Vs values for tills range from 200 to 1300 m/s as given by MASW and downhole meas-

V s30 Estimates from Borehole Data
The two datasets of the MTQ and the MELCC are used to estimate V S30 . Each selected borehole site includes information on the soil types, their relative thicknesses, as well as the RQD of bedrock. Vs values of the four main soil types (backfill, clay, sand and till) were estimated using the proposed Vs relation for the four prevalent soil layers in Montreal provided by Rosset et al. [12]. For clay and sand types, Vs profile with depth H is proposed combining various seismic measurements and borehole data. Table 3 lists the model coefficients, and Figure 6 show the data and model for sand and clay up to 35 m.

V s30 Derived from HVSR Resonance Frequency
The HVSR method is a widely used analysis technique that estimates the fundamental frequency of a site through the calculation of the ratio of the Horizontal to Vertical Fourier (H/V) Spectra obtained from ambient noise or microtremor recordings. The influence of the soft soil thickness is correlated to the fundamental frequency (f 0 ) by the general relationship f 0 = Vs/4H where Vs and H are the shear-wave velocity and thickness of the soil layer. A three-component seismometer is used to record ambient noise for twenty-minute at a sampling rate of 128 Hz, which is sufficient to obtain stable, results [18]. Each ambient noise recording was analyzed with the Geopsy tool [26] in order to identify the frequency of the maximum amplitude peak of the HVSR.
Most of the ambient noise recordings are from a 3-component sensor Tromino ® . The derived f 0 is compared to the estimated V s30 in sites at a distance less than 250 m of an HVSR site. Rosset et al. [12] developed a relation using 148 sites, mainly in Montreal, Laval and Longueuil, which is updated by 39 new sites acquired in 2022 in the neighborhood municipalities. The updated equation is of the form: The values of f 0 and V s30 at sites where both observations are available are plotted in Figure 8. The newly calculated regression equation (solid line) has a coefficient of determination R 2 of 0.8 and a standard deviation of around 50 m/s. The dashed line shows the previous equation proposed by Rosset et al. [12]. V s30 value derived from f 0 using Equation (3) is estimated for 2653 sites, as shown in the map in Figure 9. Among these sites, around 15% of them have a second frequency peak with lower H/V amplitude. and Longueuil, which is updated by 39 new sites acquired in 2022 in the neighborhood municipalities. The updated equation is of the form: The values of f0 and Vs30 at sites where both observations are available are plotted in Figure 8. The newly calculated regression equation (solid line) has a coefficient of determination R 2 of 0.8 and a standard deviation of around 50 m/s. The dashed line shows the previous equation proposed by Rosset et al. [12]. Vs30 value derived from f0 using Equation (3) is estimated for 2653 sites, as shown in the map in Figure 9. Among these sites, around 15% of them have a second frequency peak with lower H/V amplitude.

Microzonation in Terms of Vs30
A total of 229,387 sites where the Vs30 value is estimated from the different sources of data are used to produce a zonation map.
The values are interpolated in a grid of points equidistant of 0.01 degree in longitude

Microzonation in Terms of V s30
A total of 229,387 sites where the V s30 value is estimated from the different sources of data are used to produce a zonation map.
The values are interpolated in a grid of points equidistant of 0.01 degree in longitude and latitude (corresponding to a cell size of 300 m and 400 m, respectively) using the inverse distance weighted method [27] with the power parameter p fixed to 4. The selected region is a rectangle including the MMC with minimum and maximum coordinates (45.229 • N, 74.268 • W) and (45.984 • N, 73.062 • W). The interpolated V s30 are grouped into ranges of values corresponding to the NBCC site classes shown in Table 1 (Figure 10). In addition, cells of the regular grid with outcropping bedrock are directly indicated as class A (V s30 = 2082 m/s).   Figure 10. Vs30 map grouped in site classes (Table 1) using Inverse Distance Weighted interpolation procedure.

Discussion
In Figure 10, we propose the first Vs30 microzonation map for the region of Montreal using direct or indirect Vs profiles, borehole and geological data. It uses the approach developed earlier for the island of Montreal [12] by including a newly available online borehole dataset, unpublished seismic refraction data and additional sites investigated by the HVSR method. The latter is an effective non-invasive seismic survey method in terms of cost and time when sufficient data have been collected from MASW post-glacial sites to correlate the dominant frequency peak f0 with Vs30 (Equation (3)). Nevertheless, Motazedian et al. [14] indicate that in Ottawa, HVSR f0 underestimates the frequency based on the relation Vsav/4H, where Vsav is the average shear-wave velocity for postglacial sites with thicknesses greater than 50 m.
Seismic refraction profiles have three layers for which the soil type is not always specified. The data set of Vp/Vs values for similar soil types in Ontario helps us to convert the provided Vp values to Vs values for further calculation of Vs30. An improvement to refine the proposed regression equations in Table 2 could be to conduct seismic surveys at the location of multiple SR sites.
The borehole dataset is a useful complementary source of information when available but requires detailed information to correlate soil types and Vs estimates. In our case, the Vs30 uncertainty associated with this dataset is relatively high, and in several regions, Figure 10. V s30 map grouped in site classes (Table 1) using Inverse Distance Weighted interpolation procedure.
The results indicate that the Northeastern part of the MMC (Les Moulins, L'Assomption, Marguerite d'Youville and Vallée de Richelieu) is mainly in class D (V s30 between 180 and 360 m/s) and E (<180 m/s) and are mainly derived from SR and HVSR data ( Figures 5 and 9, respectively). Vaudreuil-Soulanges, in the southwest part, is similar in class D from HVSR and borehole data (Figure 7). Thérèse de Blainville, in the north, is dominated by class D sites and derived mainly from borehole data (Figure 7). Table 4 shows the distribution of site classes by data source. Site classes D and E predominate in the SR and field experiment data sets. For both of them, the selection of sites in zones where the thickest layer of soft soil was expected in Montreal Island is the first explanation for low V s30 values. Borehole sites are mainly classified in classes B and C, with 21% of the sites in class A, while the HVSR data has 23% in class D and a relatively similar number in classes B and C. Again, the latter distribution is influenced by the selection of sites where soft soils were expected.

Discussion
In Figure 10, we propose the first V s30 microzonation map for the region of Montreal using direct or indirect Vs profiles, borehole and geological data. It uses the approach developed earlier for the island of Montreal [12] by including a newly available online borehole dataset, unpublished seismic refraction data and additional sites investigated by the HVSR method. The latter is an effective non-invasive seismic survey method in terms of cost and time when sufficient data have been collected from MASW post-glacial sites to correlate the dominant frequency peak f 0 with V s30 (Equation (3)). Nevertheless, Motazedian et al. [14] indicate that in Ottawa, HVSR f 0 underestimates the frequency based on the relation V sav /4H, where V sav is the average shear-wave velocity for postglacial sites with thicknesses greater than 50 m.
Seismic refraction profiles have three layers for which the soil type is not always specified. The data set of Vp/Vs values for similar soil types in Ontario helps us to convert the provided Vp values to Vs values for further calculation of V s30 . An improvement to refine the proposed regression equations in Table 2 could be to conduct seismic surveys at the location of multiple SR sites.
The borehole dataset is a useful complementary source of information when available but requires detailed information to correlate soil types and Vs estimates. In our case, the V s30 uncertainty associated with this dataset is relatively high, and in several regions, the high density of sites could lead to a large discrepancy in the parameter, an indicator of the rapid spatial variability in the soil profiles. A preliminary analysis of the spatial correlation of V s30 derived from borehole data indicates that only short correlation distances are significant. Future spatial interpolations should account for the spatial correlations and the uncertainties associated with V s30 from each source and report best estimates, as well as their uncertainties as a function of location (Table 4).
The validation of the produced map is limited because the region of Montreal is rarely affected by earthquakes felt by many and by the lack of seismic instruments to record in several location weak motions. Nevertheless, Rosset et al. [11] have compared the distribution of site classes from three microzonation maps, including the one presented here, and the collected intensity data derived from the automatic Did You Feel It (DYFI) system operating since 2000 and from the reports collected on paper after the M5.8 1988 Saguenay earthquake. For the MMC, they used around 12,500 reports from 23 earthquake events (magnitude varying from 2.7 to 5.9) to estimate the mean felt intensity in a regular grid of cells considering the cells with more than four reports (400 cells over the 1865 cells with reports). They have shown that an increase in mean intensity of about 0.4 units with site class is observed in the MMC when moving from class A to E, suggesting that site conditions could influence the local perception of ground shaking [11]. The proposed site class's distribution is superimposed on the calculated grid of felt intensity provided by Rosset et al. [11] in Figure 11.
In general, the cells with the highest values of intensity are located in site classes E (in red) and D (in purple) of the microzonation map (see [11] for further analysis). This is not true in some areas where site classes are B (in blue) or C (in yellow) and average intensity is higher than III. In other cases, low intensity is found in areas with site classes D to E, where we could expect site amplification. In both cases, these areas should be further investigated to confirm or infirm the citizen's observations. Figure 11. Comparison of average intensity from DYFI data and the site class's microzonation using the uniform IDW procedure. The average MMI for cells with at least five reports is shown by circles with size proportional to the intensity. Values higher than III are written in the circle (from [11]). Figure 11. Comparison of average intensity from DYFI data and the site class's microzonation using the uniform IDW procedure. The average MMI for cells with at least five reports is shown by circles with size proportional to the intensity. Values higher than III are written in the circle (from [11]).
In general, the cells with the highest values of intensity are located in site classes E (in red) and D (in purple) of the microzonation map (see [11] for further analysis). This is not true in some areas where site classes are B (in blue) or C (in yellow) and average intensity is higher than III. In other cases, low intensity is found in areas with site classes D to E, where we could expect site amplification. In both cases, these areas should be further investigated to confirm or infirm the citizen's observations.

Conclusions
The next step, to account for the uncertainty associated with the different data sources, is to include a variable weight in the interpolation procedure and propose a second map showing the standard deviation. The high density of borehole data in certain regions will also be taken into account by averaging over a regular mesh the V s30 values in order to consider the large spatial variability of the V s30 value at short distances. Table 5 gives an estimate of the standard deviation to be considered in the V s30 calculation. For MASW, HRSR, and downhole sites, the estimate is fixed empirically at 30 m/s, as if the downhole method was expected to give more accurate results. For HVSR sites, the standard deviation calculated in Equation (4) is used. For SR and borehole data, the standard deviation is estimated by propagating the standard deviation given in the different equations in Sections 2 and 3 in the calculation of V s30 . The standard deviation for the sites of the grid where the rock is identified from the digital geological map [24] is fixed at 100 m/s. Information from Figure 11 will also be used to validate certain zones of the microzonation where reported intensities are not in line with expected trends, as has been done, for example, by Van Noten et al. [28] in the Anglo-Brabant Massif and the Lower Rhine Graben.
In addition, the deployment of NRCan strong motion instruments in the Montreal region as part of the National Earthquake Early Warning Network [29] will help to acquire new seismic records useful for comparing the distribution of site classes in Figure 10 and the spatial variability of their ground motion field.