Geotechnical investigations of an earthquake that triggered disastrous landslides in eastern Canada about 1020 cal BP

Abstract

Geotechnical studies are carried out to investigate the location and magnitude of a prehistoric earthquake in eastern Canada. Previous studies identified 12 sensitive clay landslides of about 1020 cal BP in the Ottawa-Gatineau region, including one of the most disastrous slope failures in eastern Canada. The landslides were hypothesized to have been triggered by an earthquake. In the current study, three of those landslides are investigated to understand their failure mechanisms. The study sites form a triangle of lateral distances of about 38, 38 and 35 km. Micro-seismic surveys, cone penetrometer tests, vane shear tests, soil sampling and laboratory testing are conducted to collect geotechnical data. Slope stability models are constructed to calculate threshold ground accelerations required to trigger the landslides. The ground accelerations are used with Ground-Motion Prediction Equations to triangulate the earthquake. The results indicate an earthquake of a minimum moment magnitude of 5.9 near Eardley, Quebec.

Introduction

Understanding ancient earthquakes helps improve building code by extending seismic data to prehistoric times (Adams 2011). It is important for long-term seismic design of critical infrastructures. Geological Survey of Canada (GSC) has been carrying out paleoseismicity research in the Ottawa-Gatineau region, which is located in the Western Quebec Seismic Zone in Canada. The seismic risk of the region ranks the third in Canadian urban areas (Lamontagne 2010). Previous studies dated a number of Champlain Sea clay landslides in the region (Aylsworth et al. 2000; Brooks 2013, 2015). The studies identified 12 landslides aged about 1020 cal BP (Fig. 1), including a Quyon River valley landslide, one of the most disastrous Champlain Sea clay failures (Brooks 2013). Those landslides were hypothesized to have been triggered by an earthquake. Brooks (2013) estimated a moment magnitude (M_w) of the earthquake to be at least 6.1. The estimate was based on Keefer (1984) and Rodriguez et al. (1999) correlations between earthquake magnitude and coverage area of landslides. The estimate has considerable room for improvement simply because the adopted correlations are based on worldwide data that do not consider the behaviour of sensitive clay failures nor do they take into account region specific ground-motion characteristics. The failure mechanisms of sensitive clay landslides are very different from other types of slope failures and deserve special attention. Also, soft sediments have ground-motion amplification effect (Bent et al. 2015) that is also not considered in the correlations. The ground-motion characteristics in eastern North America is different from that in western North America for example (Atkinson 2008). GSC therefore initiated a project to further understand the earthquake with geotechnical investigations. By understanding the failure mechanism of the landslides, the study derives threshold ground accelerations required to trigger the slope failures. The results are then used to calculate the minimum magnitude and epicenter of the earthquake. It adds knowledge to the earlier works by taking into account site specific geotechnical information and ground-motion characteristics of the region.

Fig. 1

Location map of the study area. Solid triangles mark the headscarp locations of the study landslides. Hollow triangles are the locations of other landslides of ~ 1020 cal BP (after Brooks 2015). Dotted lines are some major faults of interest in the region (after Atkinson et al. 2014, Bent et al. 2015, and Lamontagne et al. 2020)

Three out of the 12 landslides of ~ 1020 cal BP are selected for this study. They are located at Breckenridge, Low and Quyon, all in Quebec (Fig. 1). The landslides occurred in sensitive clays deposited between 13.9 and 11.5 cal ka BP within the Champlain Sea region (Dyke and Prest 1987; St.-Onge 2009). Typically, the sensitive clays are covered by a thin veneer of sandy materials on the surface and underlain by a layer of till or by bedrock where till is absent (Gadd 1986). The three sites were chosen because they are relatively far apart, fairly evenly spaced (35, 38 and 38 km from each other) and form a near optimal triangle for the purpose of this study (triangulation of the earthquake). Also, the selected landslides are among the largest Champlain Sea clay failures in eastern Canada. Thus, the effects of some geotechnical uncertainties are minimized. For example, the failure mechanism of larger landslides can be established with more certainty (Wang 2016, 2019b). Furthermore, with larger landslides, the accuracy of the reconstructed slope models, such as the local stream side slope geometry, has less effect on stability analysis than with smaller landslides. Another reason for choosing those three landslides is that the site conditions are relatively well preserved.

Wang (2016, 2019a and 2019b) reported detailed geotechnical studies of the Quyon, Low and Breckenridge landslides respectively. Due to space limit, the current paper outlines only the key information necessary for the earthquake interpretation.

Geotechnical studies of selected landslides

Breckenridge landslide

The landslide at Breckenridge is about 20 km northwest of Ottawa, Ontario (Figs. 1, 2). The Champlain Sea sediments in the area covers a deep oval-shaped bedrock basin of about 5 km by 4 km and a maximum depth of about 98 m (Crow et al. 2017). The Champlain Sea plain is fairly level with an average elevation of about 100 m above sea level (asl.). A postglacial stream network incised into the Champlain Sea deposits in the basin. Breckenridge Creek is a tributary channel to the Ottawa River (Fig. 2). There are 38 scars of sensitive clay landslides in the Breckenridge basin (Wang 2019b). Some of them are visible in Fig. 2. The ages of those landslides range from as recent as 2008 AD to older than 7000 cal BP (Brooks 2013). The study landslide of 1020 cal BP is the largest slope failure in the Breckenridge basin (referred to as the Breckenridge landslide hereinafter). The size of the landslide is about 980 m long and about 370 m wide. The depression of the scar zone is about 8 to 15 m below the Champlain Sea plain. The thickness of the sediments at the study landslide location is about 50 m to 60 m. The scar area is mostly wooded and the hummocky terrain of the debris surface is well preserved. The surrounding flat plains are mostly hayfields and pastures.

Fig. 2

Breckenridge landslide site (LiDAR image© Government of Quebec)

Geotechnical field tests at this site include Cone Penetrometer Tests (CPT) at four locations from CPT1 to CPT4, and Vane Shear Tests (VST) at three locations from VST1 to VST3 (Fig. 2). CPT1, CPT4 and VST1 are located outside the landslide zone and the rest inside. Undisturbed soil samples were taken from borehole BH1 at depth from 3 m to 18 m. The samples were tested at the Sedimentology Laboratory of the GSC for grain size and geotechnical index properties. The results indicate all fines of mostly fat clay or elastic silt (ASTM International 2017). The CPT results show that the sediments are fairly uniform across the landslide zone. The sensitivity of the material is higher (around 30 or greater) in the upper 16 m zone and lower (less than 30) below this depth. The slip surface is found to be at about 15 m depth from the Champlain Sea plain. The CPT data are calibrated with the VST results to determine the soil peak undrained shear strength (C_u). The C_u value increases linearly with depth as follows:

$$ {\mathrm{C}}_{\mathrm{u}}=15+2.66\ \mathrm{H} $$

(1)

where H is depth (m) from the pre-failure surface (100 m elevation).

Figure 3 is an interpreted cross-section along the approximate centerline of the landslide (location shown in Fig. 2). The pre-failure ground surface is interpolated from the adjacent Champlain Sea terraces (elevations from LiDAR data). The bedrock is far below the landslide depth based on the CPT results and with reference to Crow et al. (2017). Wang (2019b) presents considerable evidence and interprets the landslide as a “flake” slide (Bjerrum 1973) with a slip surface developed across the entire failure zone at the onset of the failure.

Fig. 3

Cross-section of Breckenridge landslide along the approximate centerline shown in Fig. 2 (vertical exaggeration: 2x)

A two-dimensional slope model is constructed base on the cross-section shown in Fig. 3. A computer software, Slope/W (Geo-Slope International Ltd. 2010) is used for slope stability analysis to determine a minimum ground acceleration required to trigger the landslide. Slope/W is an internationally adopted software for slope stability studies. Morgenstern-Price limit equilibrium method is used for the analysis. The limit equilibrium method is considered appropriate for the purpose of this study. For example, the method is well suited for “flake” type failures. The material parameters can be obtained with standard procedures. Although there is much room for improvement as with any model, the limit equilibrium method has high success rate with considerable case histories in the literature and in practice.

The slope model consists of one layer of soil material overlying bedrock. The surface crust is ignored for simplicity as it has negligible effect due to the great length of the slope. The stream side slope is extrapolated from the current valley side slope. The sensitivity of this assumption diminishes with the relatively long slope geometry. The slip surface is inferred from the field test results. There could be some uncertainty about the slip surface locally near the stream side, but again its effect diminishes with the relatively long slope geometry. Total stress analysis is performed with the peak undrained shear strength given in Eq. 1. A soil unit weight of 14.9 kN/m³ is assumed based on the average test results of the soil samples (Wang 2019b).

A factor of safety (FOS) of 20 is calculated for static conditions with the above slope model configuration. Obviously, it is impossible for the slope to fail statically without external loading (unless a retrogressive failure occurred which is unlikely as discussed extensively with considerable evidence in Wang 2019b).

Pseudo-static analysis is carried out to estimate the threshold ground acceleration required to trigger the landslide. A horizontal seismic coefficient is used to simulate the ground acceleration. In other words, a horizontal driving force (pseudo-static) is added to the sliding body and a factor of safety is calculated. The vertical component of the ground acceleration is ignored as it has negligible effect on the relatively level ground. The aim of the analysis is to find the ground acceleration that brings the FOS to 1.0, which is considered as a critical condition. The analysis is performed with a trial-and-error approach. A random seismic coefficient is applied to the model on the first run to calculate its corresponding factor of safety. The seismic coefficient is adjusted to another higher or lower random value if the FOS from the previous round of calculations is higher or lower than 1.0. The analysis is repeated with the new seismic load applied to calculate its corresponding factor of safety. The iteration repeats until a FOS = 1.0 is achieved. With this iteration process, a threshold ground acceleration of 0.28 g (for FOS = 1.0) is obtained that is considered as a minimum value at the Breckenridge site.

Quyon landslide

The Quyon River valley landslide (or Quyon landslide in short) is an extremely large and complex slope failure (Figs. 1 and 4). The landslide is about 35 km WNW of the Breckenridge site. It covers an area of about 14 km long and 4 km wide. The landslide was previously described as multiple separate slope failures until it was identified by Brooks (2013) as a single event that he hypothesized to have been triggered by the subject earthquake. The failure occurred along both sides of the Quyon River, a tributary of the Ottawa River. The area in and around the landslide is primarily agricultural consisting of cultivated fields, pasture and woodlands. There are two main failure zones: a northern zone and a southern zone (Fig. 4). The two failure zones are separated by a narrow area of about a kilometer long. The surface of the Champlain Sea sediments is sub-horizontal with a longitudinal gradient (along the Quyon River valley) of about 0.1° in the northern zone and about 0.5° in the southern zone. The general depression of the scar zone can be as deep as 40 m from the original Champlain Sea plain. Most of the scar zones are wooded and hummocky. The terrain of the debris surface has been well preserved in the wooded areas where horsts and grabens typical of sensitive clay landslides (Locat et al. 2011) are visible.

Fig. 4

Study area and shaded relief map of the Quyon landslide (LiDAR image© Government of Quebec. Yellow dash-dotted line marks the approximate centerline of the failure zone where longitudinal failures are analyzed. Red lines are locations where the lateral sections are analyzed)

Geotechnical investigations at the Quyon landslide site consist of micro-seismic surveys, CPT, VST, soil sampling and laboratory testing. The micro-seismic survey was conducted with a handheld device (Tromino®) to measure the sediment thickness. A total of 180 locations in and around the landslide scar were surveyed. Details of the survey locations and the results are available from Wang et al. (2015). The thickness of the sediments ranges from zero at bedrock outcrops to as deep as 68 m below the current debris surface. A bedrock contour plot indicates a bedrock valley underlying the Quyon River (Wang et al. 2015). The bedrock valley narrows down at the location between the two main landslide zones.

Drillhole tests consist of six CPTs (CPT1 to CPT6) and two VSTs (VST1 and VST2) at locations shown in Fig. 4. All CPT holes were drilled to refusal (bedrock) with depth ranging from 23 m (CPT6) to 55 m (CPT4) below the current surface. The vane shear tests were done at 5 m intervals from 5 m to 29 m depth at VST1 and from 13 m to 49 m at VST2. Undisturbed soil core samples were taken from the VST holes at depths between the vane shear tests. The samples were tested at the GSC’s Sedimentology Laboratory for gradation and geotechnical index properties. The samples are almost entirely fines (97% to 100%) according to ASTM International (2017). The materials in the southern zone are finer than that of the northern zone. The CPT data are calibrated with the VST results for peak undrained shear strength. The sediments in the southern failure zone are weaker than that in the northern zone. The peak undrained shear strengths fall within the following bounds:

$$ \mathrm{Lower}\ \mathrm{bound}:{\mathrm{C}}_{\mathrm{u}}=10+2.26\ \mathrm{H} $$

(2)

$$ \mathrm{Upper}\ \mathrm{bound}:\kern0.5em {\mathrm{C}}_{\mathrm{U}}=35+3.36\ \mathrm{H} $$

(3)

where C_u is peak undrained shear strength (kPa), and H is depth (m) from the pre-failure ground surface.

The lower bound value (Eq. 2) is more representative of the materials in the southern failure zone and the upper bound (Eq. 3) more of the northern zone.

Wang (2016) presents considerable evidence and interprets that the Quyon landslide is a “flake” slide with an overall slip surface developed at the onset of the failure. Therefore, pseudo-static limit equilibrium analyses are carried out to estimate the threshold ground acceleration required to trigger the landslide. Figure 4 shows the locations of seven sections where slope stabilities are analyzed: one longitudinal section along the approximate centerline of the landslide and six lateral or cross-sections at S1 to S6. Figure 5 is a profile of the longitudinal section and Fig. 6 is a typical lateral section (S5). The pre-failure surface is interpreted from the undisturbed Champlain Sea terraces along the shoulders of the landslide scar. The current ground (debris) surface is from the LiDAR data and is shown for reference. The bedrock surface is interpreted from the micro-seismic survey results. The slip surface is interpreted from the CPT results as well as observations along the Quyon River (Wang et al. 2015; Wang 2016).

Fig. 5

Longitudinal section of Quyon landslide (vertical exaggeration: ~16x)

Fig. 6

A typical lateral section of Quyon landslide (S5 shown in Fig. 4; vertical exaggeration: ~6x)

For the longitudinal section, the northern and southern failure zones are analyzed as two separate failures. For the lateral sections, the slope on each side of the Quyon River is modeled independently from the other. A total of 14 slope models are therefore analyzed: two for the longitudinal section and 12 for the lateral sections. One layer of soil material is assumed overlying bedrock for each slope model. The surface crust is ignored due to its negligible effect on the relatively long slope geometry. The pre-failure river side slopes are extrapolated from the current river side slope. The accuracy of this interpretation has negligible effect on the model due to the long slope geometry. The lower bound and upper bound strength properties given in Eqs. 2 and 3 are used in the analyses. The average soil unit weights of 14.9 and 17.4 kN/m³ obtained from laboratory tests are assumed for the lower and upper bound cases respectively.

Under static conditions, the lowest factor of safety of all the model cases is 5.8. Pseudo-static analysis is carried out with a trial-and-error approach as described for the Breckenridge case. A threshold horizontal ground acceleration is obtained from each model when FOS = 1.0 is achieved.

Tables 1 and 2 provide a summary of the threshold ground accelerations from the 14 slope models. As mentioned earlier, the upper bound material properties are more representative of the northern failure zone. The upper case results of the northern zone models are therefore likely the governing conditions. In other words, the highest horizontal ground acceleration of 0.27 g in Tables 1 and 2 is likely required to trigger the entire landslide. It is noted from the Tables that the threshold 0.27 g is predicted from both the longitudinal and lateral models. It is an indication that the failure could have occurred either longitudinally along the Quyon River or laterally perpendicular to the river or a combination of both. This is consistent with the complex failure phenomena of the Quyon landslide, where the LiDAR map indicates the debris deformation in directions both along and perpendicular to the river depending on location.

Table 1 Threshold horizontal ground accelerations (g) of longitudinal models of Quyon landslide

Table 2 Threshold horizontal ground accelerations (g) of lateral sections of Quyon landslide

Low landslide

The Low landslide is about 38 km north of the Breckenridge site and about 38 km northeast of the Quyon landslide (Figs. 1 and 7). The slope failed into Stag Creek, a tributary to the Gatineau River. Regionally the topography is hilly with elevations ranging from about 100 m to 250 m above sea level. Locally, the ground elevation varies from about 113 m at the creek level to about 170 m at nearby rock outcrops. Champlain Sea sediments are present in the area with a flat plain of about 153 m (asl.) surrounding some elevated sediment lobes (Fig. 7).

Fig. 7

Low landslide site (LiDAR image©Government of Quebec. White dotted line marks the perimeter of the large failure zone of 5200 cal BP. Orange dashed line is the boundary between an upper debris plateau and a lower debris plateau or headscarp of the 1020 cal BP failure. The yellow dash-dotted line is the approximate centerline where slope stability is analyzed)

There is an elongated failure zone of about 790 m long and 430 m wide (Fig. 7). The depth of the scar zone is about 10 to 15 m below the Champlain Sea plain. There are two distinct debris plateaus inside the scar zone as shown in Fig. 7. The elevation difference between the upper and lower plateaus is about 5 m. The lower zone is about 520 m long. The boundary between the two plateaus is horseshoe shaped. There is a fairly sharp slope face of about 25° (or ~ 2H:1 V) along this boundary line, which appears to be a shear failure plane.

Open land mostly hayfields and pastures surround the landslide. The area inside the landslide scar is wooded in the southeast quadrant and pasture or wetland in other areas. The debris surface appears well preserved. The terrain outside the landslide zone also appears little altered by human activities as seen from the LiDAR image.

To the east of O’Conner Road, surface water discharges to a gully and flows eastward away from the landslide zone. The drainage pattern inside the landslide zone is westward to Stag Creek. The creek is about 5 m wide and about a meter deep. The creek side slopes are mostly covered by trees. The side slopes are generally about 2.3H:1 V, with the height ranging from 10 m to 35 m. Local surface sloughing or small scale slope failures are frequent as observed from vegetation changes or freshly exposed soils. The creek is incised nearly to bedrock. Probing along the creek encountered bedrock at about 2 to 3 m depth at some locations. Large areas of horizontally stratified clays are visible along the creek bed. At some locations, the clays are exposed by some tens of meters stretching along the creek bed that appear to have not been affected by the landslide.

The geotechnical field program at the Low landslide site started with micro-seismic surveys to determine the sediment thickness in and around the landslide area. A total of 27 locations inside and outside the scar were surveyed. The interpreted sediment thicknesses range from 10 to 44 m (approximate depths to bedrock from the current ground surface). Shallower bedrocks were measured near the sediment lobes (outside the scar) indicating bedrock knolls underlying the lobes.

Undisturbed soil samples were taken at depths from 3 m to 23 m (or beyond the influence of the landslide) from boreholes BH1, BH2 and BH3 (Fig. 7). The soil samples were tested at the GSC’s Sedimentology Laboratory for grain size and geotechnical index properties. The results indicate 100% fines according to ASTM International (2017).

Cone Penetrometer Tests were carried out at CPT1 to CPT4 (locations shown in Fig. 7). Two CPTs were inside the landslide scar and two outside. All CPT holes were drilled to refusal (bedrock). Field Vane Shear Tests were carried out at VST1 to VST4 that are adjacent to the CPT holes (Fig. 7). The shear strength data collected from the VST holes are used to calibrate the CPT data. The results indicate that the peak undrained shear strength of the material increases linearly with depth as follows:

$$ {\mathrm{C}}_{\mathrm{U}}=28+1.42\ \mathrm{H} $$

(4)

where C_u is peak undrained shear strength (kPa); and H is depth (m) from the original Champlain Sea plain.

The field test data and other information unveiled that the overall landslide scar of 790 × 430 m in Fig. 7 is likely a product of two staged failures (Wang 2019a). The first failure likely occurred about 5200 cal BP to a length of about 790 m and depth about 20 m from the Champlain Sea plain (corresponding to slope failures triggered by another earthquake in the region about 5200 cal BP, Brooks 2015). The debris surface of the first failure settled at about 10 m below the Champlain Sea plain. The second failure likely took place inside the older debris zone to a length of 520 m about 1020 cal BP. The depth of the second failure is about 20 m below the older debris surface. The new debris surface settled at about 5 m below the older one. Figure 8 shows an interpreted cross-sectional profile along the approximate centerline of the landslide (location shown in Fig. 7). Details of the interpretation are discussed in Wang (2019a).

Fig. 8

Cross-section of the Low landslide (location shown in Fig. 7). The CPT locations are approximate projections to the centerline. (Vertical exaggeration: 4x)

Same as the Breckenridge and Quyon landslides, the 1020 cal BP landslide at Low is likely a “flake” type failure with an initial slip surface developed along the entire 520 m failure zone at the onset of the failure (Wang 2019a). Slope stability analysis is carried out in a way similar to that for the Breckenridge and Quyon landslides. A slope stability model is constructed based on the 1020 cal BP failure configuration shown in Fig. 8. The clay shear strength is assumed to follow the correlation of Eq. 4. The unit weight of the material is assumed to be 16.8 kN/m³, an average value of the test results (Wang 2018). A factor of safety of 6.5 is obtained under static conditions. Pseudo-static analysis yields a threshold ground acceleration of 0.19 g to trigger the 1020 cal BP failure. The result indicates that the ground shaking at Low is less intense than at Breckenridge (0.28 g) and Quyon (0.27 g). This is consistent with other observations. For example, the Breckenridge and Quyon landslides are much larger. Also, more landslides of 1020 cal BP were identified around Breckenridge and Quyon than around Low (Fig. 1).

Earthquake interpretation

The magnitude and epicenter of the earthquake are estimated with available Ground-Motion Prediction Equations (GMPE). GMPEs describe ground-motion amplitudes, e.g., Peak Ground Acceleration (PGA), as functions of magnitude, distance, and site conditions. There are three different types of GMPEs suitable for Eastern North America, i.e., simulation-based method, hybrid empirical method, and referenced empirical method (Hassani and Atkinson 2015). One from each of the three categories is used in the current study: Hassani and Atkinson (2015), Shahjouei and Pezezhk (2016), and Yenier and Atkinson (2015).

Hassani and Atkinson (2015) provides a referenced empirical GMPE for Eastern North America (referred to as HA15). The model applies an adjustment factor to a reference model of Boore et al. (2014) that is developed for other active tectonic regions. The model is applicable for soil sites.

Shahjouei and Pezezhk (2016) presents a hybrid empirical GMPE for Central and Eastern North America (referred to as SP16). The model is developed for M_w 5.0 to 8.0, distance from 2 km to 1000 km and for reference rock site conditions. It can be applied to soil sites by applying adjustment factors with procedures specified in Pezeshk et al. (2015) that further recommends procedures by Atkinson and Boore (2006).

Yenier and Atkinson (2015) proposes a simulation-based generic GMPE (referred to as YA15). It can be adjusted for use in any region by modifying relevant parameters.

The complete GMPEs and procedures appear very complex that cannot be shown here without taking a sizeable space. However, for the PGA of interest in this study, the above three GMPEs can all be illustrated in the following simple form:

$$ \mathrm{PGA}=\mathrm{f}\ \left({\mathrm{M}}_{\mathrm{W}},\mathrm{R}\right) $$

(5)

where PGA is horizontal Peak Ground Acceleration; f() is a function with adjustment factors, embedded sub-functions and conditions; M_w is moment magnitude; and R is distance from rupture (approximate epicentral distance).

Note that, in addition to the variables shown in Eq. 5, YA15 requires another variable, i.e. focal depth as an input parameter for stress effect. However, the stress effect of YA15 remains constant at focal depth equal to and greater than 10 km. According to Ma and Motazedian (2012), in most part of the West Quebec Seismic Zone where the study sites are located, the dominant focal depths are about 12 km to 22 km. The current application of YA15 therefore assumes a focal depth of greater than 10 km, hence the stress effect remains constant.

The threshold ground accelerations to trigger the landslides calculated earlier are assumed as horizontal Peak Ground Accelerations generated by the 1020 cal BP earthquake. Each of the above GMPEs is applied to find a solution for M_w and R from the PGAs at the three landslide sites.

The soil site conditions are considered in the three GMPEs with an average shear-wave velocity in the upper 30 m (V_s30) of the soils. Shear-wave velocities were measured at CPT4 at the Quyon site (Fig. 4). Details are available from Wang et al. (2015). The V_s30 is calculated to be 199 m/s at the Quyon site. Similar shear-wave velocity profiles were reported for the Breckenridge site by Crow et al. (2017). Shear-wave velocities are not measured at the Low site. The average V_s30 = 199 m/s from the Quyon site is used for all the three sites for simplicity when applying the GMPEs (sensitivity check within a reasonable range of V_s30 indicated negligible effect on the results).

A solution from each GMPE can be obtained with a trial-and-error approach using any readily available spreadsheet software. A chart as shown in Fig. 9 is used in the current exercise to visualize each GMPE. A curve is plotted on the chart to show the relationship between M_w and PGA for any given distance, R. The procedure starts from specifying a random magnitude, e.g. M_w = 5.0, and plotting it as a horizontal line on the chart (Fig. 9). The next step is to look for a distance, R, that yields a M_w ~ PGA curve intersecting point (M_w, PGA), e.g., (M_w 5.0, 0.28 g), (M_w 5.0, 0.19 g), (M_w 5.0, 0.27 g) for Breckenridge, Low and Quyon sites respectively. This is achieved by curve fitting (trial-and-error). The resulting distance, R, for each curve (corresponding to each PGA or site) represents the approximate epicentral distances to that site. A location plan as shown in Fig. 10 is then used to locate the epicenter. On this chart, a circle is plotted with each landslide site as a center (origin) and the distance, R, determined above as radius. If the three circles do not cross at a single point, the specified M_w is deemed as impossible. For example, if the three circles (arcs) cross each other at three different locations, each crossing location is an epicenter that satisfies the conditions at only two landslide sites, but not the other. The assumed magnitude therefore needs to be adjusted and another round of calculations carried out. The iteration continues until the three arcs cross at a single point where the assumed magnitude satisfies the conditions at all the three landslide sites. The crossing point closest to the three landslides is considered as the epicenter and resulting magnitude is considered as the minimum magnitude of the earthquake. Figures 9 and 10 show the end results of the iterations (epicentral distances rounded to the nearest km). Table 3 summarizes the results. The magnitude is estimated to be M_w 5.3, 5.9, and 5.9 with HA15, SP16 and YA15 respectively. It happens that the estimated epicenter is located inside the triangle formed by the three landslide sites as shown in Fig. 10 (note that the procedure does not restrict the epicenter to be inside the triangle). All three GMPEs predict an epicenter in the vicinity of each other and located near Eardley, Quebec (Figs. 1 and 10) where a major fault (Eardley Fault) exists. It must be noted that it is not the intention of this paper to suggest that the deduced earthquake was associated with the particular fault. The cause of the earthquake and the region’s geological structures are complex issues that are beyond the scope of this paper.

Fig. 9

Relationship between M_w and PGA for any given distance from rupture with three Ground-Motion Prediction Equations: (a) HA15, (b) SP16, and (c) YA15. (Magnitudes are rounded to the first decimal point. Epicentral distances are rounded to the nearest km)

Fig. 10

Earthquake triangulation with three GMPEs: (a) HA15, (b) SP16, and (c) YA15 (Magnitudes are rounded to the first decimal point. Epicentral distances are rounded to the nearest km. Dotted lines mark the nearby Eardley Fault shown in Fig. 1)

Table 3 Summary of estimated earthquake magnitude and location

Discussions

The main focus of this paper is using a geotechnical approach to decode the information about the earthquake preserved with the landslides. With the resulting ground accelerations from the landslide studies, two of the GMPEs (SP16 and YA15) predict a moment magnitude of M_w 5.9 and the other (HA15) predicts a magnitude of M_w 5.3. While it is beyond the scope of the current paper to discuss the Ground-Motion Prediction Equations, Yenier and Atkinson (2015) does point out that there is “an important limitation of the referenced empirical approach” (as used by Hassani and Atkinson 2015 or HA15), i.e. “the available ground-motion data in the target region may not sufficiently represent all important regional characteristics.” Yenier and Atkinson (2015) “takes advantage of key concepts from both the hybrid empirical and referenced empirical approaches to develop a robust simulation-based generic GMPE.” For this reason, and also note that YA15 and HA15 share the same co-author, the above result of M_w 5.9 is likely more reliable than M_w 5.3. Other evidences as discussed below also support this argument.

A recent Val-des-Bois (Quebec) earthquake of June 23, 2010 helps appreciate the above estimates. The epicenter of the 2010 earthquake is about 40 km ENE of the Low landslide and about 40 km north of the Ottawa River. Its moment magnitude is M_w 5.2 (Ma and Motazedian 2012). The earthquake is known to have triggered two landslides in the Champlain Sea sediments, one at Notre-Dame-de-la-Salette about 13 km southwest from its epicenter and another at Mulgrave et Derry about 18 km south of its epicenter (Perret et al. 2013). Note that among the 12 landslides affected by the 1020 cal BP earthquake, the farthest is at Beacon Hill in east Ottawa (Fig. 1), which is about 44 km away from the currently estimated epicenter. The number of landslides of 1020 cal BP and their affected area far exceed those affected by the recent 2010 M_w 5.2 earthquake. It is therefore reasonable to believe that the estimate of M_w 5.9 by both SP16 and YA15 is more likely than M_w 5.3 by HA15. Among the 12 landslides of 1020 cal BP, 9 of them are distributed along the line between Quyon and Breckenridge (Fig. 1). Location wise, the interpreted epicentre agrees well with the concentration of the landslides.

Brooks (2013) calculated a circular area of 3300 km² enclosing the 1020 cal BP landslides. With the magnitude to landslide-area correlations by Keefer (1984) and Rodríguez et al. (1999), Brooks estimated a minimum magnitude of M_w 6.1 for the 1020 cal BP earthquake. The estimate of M_w 5.9 with SP16 and YA15, even though slightly smaller, is fairly close to Brooks estimate. The current study is perhaps more precise given that site specific geotechnical conditions and the regional ground-motion characteristics have been taken into account.

The Saguenay (Quebec) 1988 earthquake discussed in Rodríguez et al. (1999) is worth noting. In terms of the geographical location of the seismic event and the sensitive clay landslides triggered, the Saguenay earthquake is the only event among the 76 cases in Rodríguez et al. (1999) that is most comparable to the present study case. The moment magnitude of the Saguenay earthquake is M_w 5.8. The landslides found in an area of 45,000 km² are classified by Rodríguez et al. (1999) as both “disrupted” and “coherent” slides. The magnitude is smaller than the prediction from the area of affected landslides. It is also smaller than predicted by the correlations with epicentral distance Rodríguez et al. (1999). When compared with Brooks M_w 6.1 estimate for the 1020 cal BP event with the same correlations, the smaller magnitude M_w 5.9 of the current prediction agrees with the Saguenay data.

Although Keefer (1984) and Rodríguez et al. (1999) correlations cannot be used to estimate epicentral locations, their worldwide data help increase the confidence of the current interpretation. Rodriguez et al. (1999) presents correlations between moment magnitude and epicentral distance to the farthest landslides of three different categories: “disrupted”, “coherent” and “spread and flow” slides. Note again that the Beacon Hill landslide of 1020 cal BP is about 44 km ESE of the estimated epicenter. For this distance, the correlations of Rodríguez et al. (1999) indicate a moment magnitude of M_w 5.6, 6.1 and 6.3 for the three categories of landslides. The current result of M_w 5.9 by SP16 and YA15 falls within this range. It is another indication that the estimated magnitude and epicenter are fairly reasonable.

Conclusions

Geotechnical data unveiled that the ground accelerations at the Breckenridge, Low and Quyon sites are at least 0.28, 0.19 and 0.27 g respectively when the earthquake occurred. The previously hypothesized earthquake associated with the landslides is interpreted with available Ground-Motion Prediction Equations using the derived ground accelerations. The results indicate a minimum moment magnitude of M_w 5.9 and the epicentre near Eardley. The findings are supported by a higher concentration of landslides in the vicinity of the predicted epicentre. The results are also consistent with observations from some recent earthquake events in the region and other world wide data. With geotechnical data, slope stability and the region-specific ground-motion characteristics considered, this study brings the understanding of the 1020 cal BP earthquake to another level.