Development of a methodology for predicting landslide hazards at a regional scale

Landslide risk analysis is a common geotechnical evaluation and it aims to protect life and infrastructure. In the case of sensitive clay zones, landslides can affect large areas and are difficult to predict. Here we propose a methodology to determine the landslide hazard across a large territory, and we apply our approach to the Saint-Jean-Vianney area, Quebec, Canada. The initial step consists of creating a 3D model of the surficial deposits of the target area. After creating a chart of the material electrical resistivity adapted for eastern Canada, we applied electric induction to interpret the regional soil. We transposed parameter values obtained from the laboratory to a larger scale, that is to a regional slope using the results of a back analysis undertaken earlier, on a smaller slide within the same area. The regional 3D model of deposits is then used to develop a zonation map of slopes that are at risk and their respective constraint areas with the study region. This approach allowed us to target specific areas where a more precise stability analysis would be required. Our methodology offers an effective tool for stability analysis in territories characterized by the presence of sensitive clays.


Introduction
Landslides affect numerous areas around the globe, damage infrastructures, and cause the loss of human lifes. This phenomenon occurs in various soils, including sands and clays. When they occur in clay soils, and depending on the sensitivity of the clay layers, landslides can become unpredictable and affect large areas. Hazard assessments thus represent vital tools in identifying and minimizing landslide risk. These assessments are based on the stratigraphy and the geometry of the different soil layers, soil types, and the assessment of geotechnical parameters, such as the shear strength of the soils ).
An initial series of key factors of landslide risk analysis involves slope geometry and terrain stratigraphy, which combine several elements that control landslide behavior. These elements include various aspects of ground topography, such as the inclination and the orientation of the slope surface and bedding planes, the presence and depth of incision by surface streams, and the proximity of the bedrock (Geertsema and L'Heureux 2014). Combining a sufficiently precise topography with a reliable description of the different soil layers permits creating an accurate 3D model of the zone in question. The accuracy of this 3D geological model is critical for an appropriate assessment of landslide and other geotechnical risks at a regional scale. Topographic information is obtained from existing maps, from geostatistical interpolation based on site surveys, and from LiDAR (light detection and ranging) surveys. Stratigraphic data is acquired through intrusive methods, such as drilling, or non-intrusive methods, such as geophysical surveys. Drilling provides core samples, which are used to identify the nature of the different soil units and serve for further analyses. Surface geophysical methods identify the stratigraphy; these methods include seismic refraction, electrical resistivity tomography (ERT), and transient electromagnetic (TEM) methods (Mussett and Khan 2000).
Of these abovementioned methods, core sampling is the most accurate but also the most expensive. Given that the required 3D model must often cover a large territory, a more suitable approach involves a large quantity of quick and inexpensive surveys across the area. When choosing the investigation methods, it is important to select those best suited to the specific area under study. In the case of a zone characterized by a great thickness of sensitive clay, a powerful detonation-required for seismic refraction-could trigger a slide and is therefore not appropriate. Electrical resistivity can differentiate clay from sand, gravel, and bedrock; this approach is useful in areas where sensitive clay layers must be differentiated from other geological units (Mussett and Khan 2000). This technique also presents the advantage of being able to detect variations in porewater salinity of the water, which affects the structure of clay. Finally, TEM is often favored because of its practicality. It is easy to transport, its installation is fast, the duration of the survey is relatively short and it can be used in all kinds of conditions. Nonetheless, interpreting these geophysical surveys can be difficult and requires adapting the interpretations to the specific target area.
A second series of key factors of landslide hazard assessment includes landslide type and the potential extent of the slide. Landslide type depends on the soil type. Multiple sophisticated methods have been developed to evaluate the possibility of landslides in clayed soils and other granular materials. These approaches include analytical and numerical stability analyses using SRF (strength reduction factor) for estimating a single, commonly spherical, failure surface landslide (Richer et al. 2020).
Extensive areas in eastern Canada are covered by sensitive clays. These clays are fine-grained plastic and cohesive soils that may or may not contain clay minerals (Holtz and Kovacs 1981). Sensitivity is an important parameter and is defined as the ratio of the undisturbed soil strength (S u ) to the remolded shear strength (S r ) . The "sensitive" qualifier is assigned to a soil having a sensitivity (S t ) greater than 1. Sensitive clays were formed in a post-glaciation context. In North America, the retreat of the Laurentide ice sheet, which began between 18 and 13 ka BP depending on the region (Dyke and Prest 1987), was followed by a marine invasion along much of the isostatically depressed coastline. The clay deposited on the coastal seafloor contained an important quantity of polyvalent cations that flocculated to form a relatively strong structure. Following the isostatically forced regression, the clay was exposed to leaching by fresh meteoric water. This leaching increased the interparticle repulsive forces.
According to the classification of Varnes (1978), and its modification by Hungr et al. (2014), sensitive clays can produce a progressive spread or a flow landslide, also called a flowslide. Both categories of landslides are highly retrogressive (Mitchell and Klugman 1979;Marko et al. 2010;Demers et al. 2014). Researchers claim that sensitive clays with a liquidity index (IL) of more than 1.2 are susceptible to retrogressive landslides (Leroueil et al. , 1996Locat and Demers 1988), although it should be emphasized that a liquidity index (IL) of more than 1.2 is only possible when the remolded shear strength (Cur) of the clay is less than 1 kPa (Tavenas 1984;Thakur et al. 2014). According to Leroueil et al. (1996), the extent of retrogressive movements of slopes and liquidity index (IL) are positively correlated. Additionally, it is evident from several slopes in eastern Canada that no retrogressive landslides occurred even when the liquidity index was above 1.2 . However, liquidity index (IL) along with other interrelated parameters such as remoulded shear strength (Cur) and quickness (Q) can be used to describe the flow behavior of landslides (Strand et al. 2017).
Given the complexity of landslides in sensitive clays, determining the landslide constraint zones is commonly based on geometric criteria. The primary methods for zoning retrogression distance include the 1:15 method (Haugen, Tveit, and Heyerdahl 2017), the "natural hazards-infrastructure for floods and slides" method (Natural hazards-infrastructure for floods and slides 2016), and the Norwegian Geotechnical Institute (NGI) method (Haugen, Tveit, and Heyerdahl 2017;Turmel et al. 2018). These approaches establish the limit of the retrogression distance at X times the height of the slope, with X varying from 2 to 15 depending on the method and the soil properties. The distance is calculated from the toe of the slope or from the base of the critical failure surface. In Quebec, a commonly used method is based on a statistical analysis of previous retrogressive landslides in the same area (Turmel et al. 2018). All these methods produce conservative results but do not cover worst case scenarios, which are probabilistically quite rare.
The third series of key factors is the accurate estimation of geotechnical parameters in the laboratory and estimating these parameters at the actual in situ scale.
The geotechnical parameters of the different soil layers produce large variations in behavior between a sample of collected material and that of an entire slope. A slope is commonly assumed to have a homogeneous composition; however, it is impossible to confirm the perfect homogeneity of a slope even in cases where a large number of samples are available. Therefore, the scale effect has an important influence on the geotechnical parameters. One of the best approaches for determining the effective value of a soil parameter is through the back analysis of past landslides (Saeidi et al. 2016). A back analysis can estimate the effective geotechnical parameters of the soil layer, where landslide scars are visible in the study area, and a reliable pre-and post-slide terrain geometry can be determined.
In this paper, we propose a methodology for predicting a regional-scale landslide occurrence and zoning a territory covered by sensitive clay slopes. This approach includes new developments in regard to two key aspects of predicting landslide occurrence: reliability of the 3D model and the assessment of geotechnical parameters. First, we obtain a reliable 3D geological model using geophysical methods and an adapted chart of the material electrical resistivity developed for a sensitive clay-covered territory. This chart is used to establish the stratigraphy of the territory and to create a 3D model of the soil deposits. We then obtain more precise values of the soil parameters through a back analysis of an earlier slide in the study area combined with laboratory data from soil samples collected from the study slide. The resulting geological model and geotechnical data are applied to the mapping of the atrisk area, and we use the method proposed by the Quebec (Canada) Ministry of Transport, the ministry responsible for evaluating terrain stability (MTMDET 2017;Turmel et al. 2018). Slopes at risk can then be identified. We identify the different stages of a regional landslide risk analysis, describe a step-by-step procedure for this analysis, and highlight relationships between the various stages.
In the following sections, we present the selected applicable methods for acquiring data and provide a general methodology for assessing regional landslides. We then describe the specific study area to which we apply our approach and present and discuss the results and implications. Figure 1 provides an overview of our methodology for determining a regional landslide risk. We detail the core procedures (blue rectangles in the figure) in the following sections.

Creating a database
First, it is important to conduct an inventory of all relevant data available for the study area. This essential information includes regional soil types, their stratigraphy, and their physical, mechanical, and geotechnical characteristics. Borehole data, geotechnical surveys, and research reports are relevant data sources.
With the existing data in hand, it is easier to identify missing information. Geophysical methods, such as the transient electromagnetic (TEM) method proposed here, can be applied to fill these gaps. TEM is an easy-to-perform non-intrusive electromagnetic induction technique. Control points must be identified and are generally borehole logs, which provide a representative stratigraphy of the area. TEM surveys are first conducted near control points, and the results are used to calibrate survey results obtained elsewhere across the study area.

Developing the 3D model
Prior to undertaking the TEM surveys, we produce an initial regional stratigraphic model using the known data points. This stratigraphic model can then be used to calibrate the collected TEM data and thus obtain a more accurate estimate of the electrical resistivity of each layer.
Once the TEM surveys are completed, we compare the resistivity values obtained through the stratigraphy of control points selected from the known data points. The objective of this step is to associate a range of resistivity values to each layer described in the borehole logs. The Palacky chart (Palacky 1987) can also serve as a reference to identify aberrant values (Fig. 2). This chart is adapted to the specific study region to identify soil types and is then used to interpret the stratigraphy recorded by all TEM surveys.
With all this data in hand, we can then overlay a grid over the study area to establish section lines. The number of created sections depends on the amount of data and the desired precision. For each section, the stratigraphies obtained from borehole logs and TEM surveys are displayed. Geological knowledge related to the deposition of each soil layer is applied to interpreting the stratigraphy between data points obtained from borehole logs and TEM surveys. At every point of intersection between two sections, both sections must display an identical interpretation of the stratigraphy. Once all sections are properly interpreted, the data set can be imported into 3D modeling software to create the 3D model.

Assessing the values of soil geotechnical parameters
Laboratory tests should be conducted on soil samples collected from each of the main soil types within the study area. The objective is to estimate the shear strength parameters (friction angle and cohesion) of each soil type and, for clays and their physical parameters like the plasticity and liquidity indices. The obtained results serve as the basis for the back analysis described below. The developed 3D model combined with site visits can identify the scars of past landslides within the study area. Pre-slide geometry is determined through the use of sequence aerial photography, LiDAR surveys (Gouvernement du Québec 2021), previous research, and the differential localization of known structures both pre-and post-slide (Saeidi et al. 2016). These marker points are then entered into a conceptual model of the slide area. Analytical or numerical methods can then be used, and realistic values are assigned to all soil materials present on the slope. The aim is to create the original slide conditions and obtain a failure surface that is similar to the observed failure surface with a safety factor of 1. The failure conditions attained in the model should produce estimates for parameters such as the friction angle and the cohesion at the time of failure. Once the conditions of the actual failure are determined, the obtained values for the slope geotechnical parameters are relatively realistic.

Developing the landslide zonation map
A landslide zoning map can be constructed to delimit the areas subject to retrogressive landslide hazards. A mapping method has been proposed for eastern Canada by MTMDET (2017) and Turmel et al. (2018); however for other locations, as Norway, various methods can be used, e.g., the 1:15 and the NIFS approaches (Haugen et al. 2017;Karlsrud et al. 1985;NIFS 2016).
Topography is used to identify all slopes at risk due to their geometry (height ˃5 m and an angle of inclination ˃14°). These slopes are then characterized further by their most representative material and the presence or absence of erosion. A realistic assumption is that landslides that have occurred over time in a given area result from specific conditions, and areas sharing similar characteristics should produce landslides having similar retrogression distances. Thus, a statistical analysis of previous landslides in the area can be used to determine the potential distance of retrogression for each slope. This approach can determine constraint areas as potential danger zones at the slope top and the bottom.

Conducting a site stability analysis
Detailed stability analysis can be undertaken for specific needs, i.e., types of infrastructure, on sites located within a constrained area. More advanced numerical modelling of flowslides and spread landslides can also be undertaken.

Study area
We selected the Saint-Jean-Vianney (SJV) area in the Saguenay region of Quebec, Canada, to test our methodology. This area is well known for its very thick deposits of sensitive clay and contains the scars of at least two major landslides, which occurred in 1663 and 1971. The 1663 landslide affected approximately 22 km 2 (Lasalle and Chagnon 1968). This event is at the origin of the designation of les terres rompues ("broken lands") for this area along the Saguenay River (Bouchard 1991). The 1971 landslide was much smaller. On May 4, 1971, during a particularly rainy spring, part of the village of Saint-Jean-Vianney was swept away by this landslide. A total of 42 houses were destroyed, and 31 people died during the event.
A geotechnical study carried out following the 1971 event (LaRochelle 1974) characterized the main parameters of the landslide. However, given the specific context of this landslide inside a larger and older slide, it remains difficult to assess the possibility of a similar event reoccurring. To commemorate the 1971 tragedy, regional authorities are planning installations for educational, geotourism, and historical purposes in the slide area (Ville de Saguenay 2007). The security of the site and its users requires that a geotechnical study be conducted over the larger area surrounding the scar of the 1971 event.

Developing the 3D model of the SJV area Creating the database for the SJV area
A hydrogeological database has been developed for the Saguenay-Lac-Saint-Jean region as part of a groundwater data acquisition program-Programme d'acquisition de connaissance sur les eaux souterraines; CERM-PACES (2013). This regional database includes data from all regional boreholes and wells, other relevant information, and the interpretation of stratigraphic sections (Fig. 3A). The stratigraphic Sects. 511 and 514 from this CERM-PACES database indicate that a clay layer, averaging 30 m thick, covers the bedrock. An intercalated discontinuous till layer, 1 to 5 m thick, is also present.
Given that the objective of the PACES project was related to groundwater resources, the thick clay deposits of the SJV area were not fully described in the initial CERM-PACES (2013) project, apart from the abovementioned sections.
For topographic data, we relied on a recent LiDAR (light detection and ranging) survey, taken in 2017, covering the SJV area. Since the occurrence of the 1971 landslide and the stabilization works conducted in the later years (LaRochelle 1974), no major activities have occurred in the area. The sole exception would be the use of the sector by all-terrain vehicle enthusiasts. Therefore, no major changes in topography were expected since the LiDAR survey was caried out (Gouvernement du Québec 2021). We undertook an initial interpretation of the regional stratigraphy with this LiDAR data. We then selected points for which the stratigraphy was known and added this new data to the database.

Acquiring and interpreting the geophysical data across the SJV area
In 2019, we conducted a field campaign using TEM to improve the data coverage within the study area. We added 75 new data points (blue dots in Fig. 3B), resulting in 18 pseudo-sections of material electrical resistivity. The control points were selected after their location near the target area for the TEM survey and their accessibility in order to optimize the time spent in the field. In this section firstly we bring the methodology of measuring and interpretation of geophysical data following by attribution TEM range of each sediment type in the region. Material and method TEM field setup A time-varying current is injected into a transmitter loop (Tx) of electrical wire deployed on the ground. This loop induces a primary electromagnetic field in the subsurface geological materials. The interaction between this field and the subsoil produces a secondary electromagnetic field that contains information on the underground electrical properties. This secondary field is detected by a receiving loop (Rx) also deployed on the ground, inside the Tx loop (in-loop configuration). The induced voltage is measured by the Rx loop, which is connected to a receptor in the TEM instrument that measures the rate of decay of the electromagnetic field. This decay rate will then be reversed in electrical resistivity (Nabighian and Macnae 1991;Fitterman and Labson 2005). The NT-32 transmitter is the TEM instrument used for this study. The NanoTEM equipment consisted of a portable battery-operated transmitter-receiver (TX-RX) console, a 32II multifunction GDP-Receiver and a high-speed sampling card with a fixed gain stage of Å ~ 10 (Mac-Innes and Raymond 2001). The TEM transmitter (TX) and receiver (RX) have a square-sized loop configuration of 20 m X 20 m and 5 m X 5 m, respectively. The filter frequency was set at 32 Hz with a 50% duty cycle, and the pulse current in the Tx loop was set at 3A with a turn-off time at 1.5 μs and a damping resistor set at 250 Ω. The data were stacked in 8 blocks of 4096 cycles each during the acquisition, giving a total of 32,768 stacks. This resulted in a noise level of approximately 10 μV/A (Fig. 4). TEM data inversion Through different software packages, the main objective is to invert the data to deduce the subsurface resistivity distribution. Converting the TEM signals into an electrical resistivity model of the ground involves three steps: (1) averaging raw data using TEMAVG Zonge software (MacInnes and Raymond 2001;2005). To significantly decrease the noise, inconsistent data points need to be filtered and deleted in this step (Fig. 4). The magnetic field intensity decreases with depth and the noise level (distortion) can increase considerably. Consequently, the decay rate of the signal sometimes becomes discontinuous at higher depths, and the results can no longer be reliable because the amplitude of the measured data decreases over time, with time playing the role of pseudodepth. To prevent the distortion of the signal, it is necessary to manually remove these inconsistent data points (i.e., red crosses (raw data) on Fig. 4) to retain only highquality data (i.e., black crosses (raw data) and black line (inversion best fit) on Fig. 4) when processing data with TEMAVG. In the first step, the results are represented by a variation of the apparent resistivity (Fig. 4); (2) After importing the measured sounding curve from TEMAVG, generating a 1D inversion model of the transient EM sounding curves through STEMINV software (MacInnes and Raymond 2005). Based on the iterative Occam inversion scheme (Constable et al. 1987). STEMINH is used to produce a consistent 1D smoothed inversion model of electrical resistivity versus depth; (3) Finally, the final phase of the data inversion allows us to build a 2D model with MODSECT software by using the 1D resistivity model acquired with STEMINV (MacInnes and Raymond 2005). MODSECT provides a 2D view to visualize the geometry of the geoelectrical structure of each TEM line by interpolating vertical columns (i.e., each TEM station) with Catmul-Rom splines.

Attribution a resistivity values for sediments in the region
As mentioned, the control points represent log data where the stratigraphy is known. During our field work, we conducted TEM surveys near these points. The control points were selected after their location near the target area for the TEM survey and their accessibility in order to optimize the time spent in the field. We then compared the electrical resistivity values for each lithology, using the Palacky (1987) chart (Fig. 2) as a reference, with the observed stratigraphy at the control points (Table 1). " The observed resistivity values for sand and gravel (Table 1) are much lower than the Palacky (1987) values. These large differences underline the importance of adapting the chart to the specific territory under study. We therefore proposed new ranges of electrical resistivity for each material. It must be noted that fieldwork occurred mostly in the spring when the regional soils are highly saturated.
In the Palacky (1987) chart, the resistivity ranges of the materials sometimes overlap; however, we preferred to differentiate each range of values and avoid any overlap. A given value fits only one material; this approach therefore facilitates interpreting the data in the subsequent steps. The minimum value for a given material was established as the average of the maximum of the lower category and the minimum of the following category. To separate clay from sand, for example, we obtained a value of 50 by rounding roughly the average of the values of 35.37 (minimum value for sand) and 72.5 (maximum value for clay). Because of a lack of values for rock, the Palacky (1987) chart set this range at a value of ≥ 1000.
Using this new chart and the pseudo-sections, we converted all TEM surveys into stratigraphy logs, and the results were entered into the database.
Interpreting the soil stratigraphy within the SJV area We interpreted the stratigraphy between data points using ArcGIS (ESRI 2015) and the Arc Hydro Groundwater tools (Aquaveo 2019).
We constructed 20 stratigraphic sections (Fig. 5) in the study area after accounting for the spatial distribution of the data points. Two sections rely exclusively on borehole logs and serve as references for the other 18 sections based on the TEM pseudo-section data. For each section, topography and the PACES-interpreted surface deposits (Rouleau and Daigneault 2013) are displayed, as well as all streams, intersections with other sections, TEM and borehole logs, and other structural elements. We used all this information to interpret each section and ensured that the intersections between sections were concordant.
During this process, we reinterpreted some TEM logs to provide a better fit to the surrounding data. As well, we found large differences in precision between the borehole and the TEM logs. Borehole data scans identify different layers of thickness at a meter-scale resolution and less, whereas TEM readings produce average values that cannot be interpreted at less than a meter scale. Given that the SJV area has been affected by at least two major landslides, these events resulted in many alterations of the thin clay and sand layers across the study area that the TEM surveys could not differentiate; we therefore had to simplify the description of some sections.
A limited number of borehole logs describe the till layer (12 of 218 boreholes), and it is difficult to differentiate till from sand or gravel on the basis of their respective electrical resistivities. Consequently, the till deposit has been integrated into the granular layers. Moreover, the low precision of the TEM survey at less than a meter-scale resolution and the complex stratigraphy of the SJV area restricted our interpretation to four layers. Ordered from the surface downward, the first is a granular alluvial sediment-dominated layer (Granular 2). The second layer is the clay-silt dominant surface layer (Clay), and the third is a granular glaciomarine sediment-dominated layer (Granular 1). The final layer is the bedrock. This simplified stratigraphy is considered sufficient for geotechnical purposes.
Creating the 3D model of the SJV area To export data from ArcGIS to Leapfrog, we used some of the methodology developed by Chesnaux et al. (2011) for building a geodatabase to map the hydrogeological features and produce a 3D modeling of groundwater systems. This approach consists of creating virtual boreholes at 50 m intervals within the previously developed sections on the basis of the geological interpretation of these sections. The augmented data set was imported into Leapfrog Geo software (Seequent 2020). The resulting 3D model consists of, from the bottom to top, the bedrock, the Granular 1 layer, the clay layer, and then the Granular 2 layer (Fig. 6).

Assessing the values of geotechnical parameters across the SJV area
During fieldwork, we collected clay cores at three locations (Tables 1, 2 and 3 on Fig. 7). We collected six Shelby tube cores, two at each location. The Shelby tubes were 700 mm long, having a 75 mm external diameter and 73 mm internal diameter. We then produced 31 samples, each approximately 10 cm long, from these Shelby cores. We used nine in the laboratory tests (Table 2), and the 22 samples were sealed in paraffin. We used six of nine samples to estimate the plastic and the liquid limits according to BNQ-2501-090; BNQ (2019) ( Table 3). The results are reliable, and average values have been computed for each location. We also carried out a sedimentation test (BNQ-2501-025; BNQ (2013) on TM2-06 and TM3-02; the results are presented in Table 4. Both samples used for the sedimentation test show near equal quantities of clay and silt (Table 4). The difference between these two samples relates to the presence of different numbers of thin silt beds. In general, the soil in the area can be defined as "CL" according to the Unified Soil Classification System (Holtz and Kovacs 1981), i.e., inorganic clays of low to medium plasticity.
Regarding geotechnical parameters of clay, nine samples were used for triaxial compression tests according to the ASTM standard (ASTM D2850-15 2015). The results presented in Table 5 include the maximum principal stress at failure (σ 1failure ), the minimum principal stress at failure (σ 3failure ), and the undrained shear strength (c u ). These values are in the range of eastern sensitive clay of Canada.
Finally, the obtained values of physical parameters and the undrained shear strength (c u ), confirm that this soil is the same type that Lefebvre (1981) have studied for determining the geotechnical parameters in the consolidated-undrained condition. They determined that the average cohesion of the clay is 8 kPa, and the average internal friction angle is 32°. Since the behavior of the soil in this occurred landslide are consolidated-undrained, these results will be used as starting point for doing back analysis for the determination of in-situ geotechnical parameters.

Back analysis of a past landslide in the SJV area
For the back-analysis phase, we selected a small landslide along the Aux Vases River to determine soil geotechnical parameters at a large scale. This small slip, which occurred between 1971 and 1975, was identified by analyzing available aerial photos covering the last 50 years.
Only two aerial photos are available for this location during this period, one in 1971 and one in 1975; thus, this landslide occurred after the 1971 SJV event and before the 1975 air photo survey.
Photo A (Fig. 8A) shows the slope condition in 2015. Photo B (Fig. 8B) shows the pre-landslide slope. The scar of an initial slip is delimited by a yellow line on Photo B. The apex of the initial slip is 40 m from the river, whereas that of the final slip is 120 m distant from the river. The elevation at the top and bottom of the slope of the initial slip is approximately 40 m and 7 m, respectively. The 40 m estimate is extrapolated from the heights of the  (Fig. 9).
We introduced the 1971 slope geometry (Fig. 8) and the present-day slope profile into the Slide software. An analytical slope stability approach (Richer et al. 2020) is used in a back analysis (Rocscience 2018). The density (γ) of the clay soil was estimated at 18.6 kN/m 3 from laboratory measurement. To estimate the shear strength parameter of the soil, several iterations of analytical slope stability modeling were carried out using various input values of soil cohesion. The initial parameter for iteration were based on the results of triaxial consolidateundrained tests conducted on in sensitive clay from Quebec, with values of 8 kPa for cohesion and 32° for friction angle (Lefebvre 1981). To determine the shear strength of the soil, we varied the values of one parameter at a time, either the cohesion value or the friction angle, as back analysis is not possible with two unknown variables. Considering the usual low cohesion values of sensitive clay in Eastern Canada, the laboratory estimated value (8 kPa) was considered as the known value.      Figure 10 illustrates the model geometry used in the Slide software and presents all possible failure surface that the safety factor is equal or less than 1. These results have been obtained after conducting a lot of iterative analysis using several values of the friction angle parameter, resulting in the nearest approximation of the failure surface for the interpreted topography. We found that the best conceptual analysis method for representing the slip profile is the "Corps of Engineers #2" (U.S. Corps of Engineers 2003). For all failure geometries that we considered, a friction angle of 35° is required to reach a safety factor of 1, a friction angle value that is near to laboratory results. Figure 10 shows that the slips profile having the minimum safety factor (SF ≤ 1) is deeper than that interpreted according to the profile of the observed present slope (purple line). The resulting shape of the slip profile is very similar to the present-day slope profile. Therefore, a cohesion of 8 kPa and a friction angle of 35° using the "Corps of Engineers #2" conceptual analysis produced the best profile of the slip.

Developing a landslide zonation map for the SJV area
The Quebec Ministry of Transport has mapped most of Quebec's territory to identify areas at risk of landslides. These zones are delimited by taking into account the characteristics of the soil and conditions related to the land, such as the soil type, the embankment slope, and the history of landslides in the zone. The Saint-Jean-Vianney area is unfortunately not mapped.
The City of Saguenay has, however, produced its own stress maps to identify zones having embankments with a landslide risk (Fig. 11). For these maps, embankment size was determined using the criteria of a height ≥ 5 m and a slope ≥ 14° (Ville de Saguenay 2007).
Using the historical aerial photos and the existing LiDAR topography, we undertook an analysis of landslides that have occurred along the Aux Vases and the  The modeled geometry of the slope in 1971 (in meters) before the small landslide Petit-Bras rivers (#1 and #2 in Fig. 11). From this analysis, we could parameterize these landslides and, in particular, estimate retrogression distance. We studied 30 landslides, including the 1971 SJV landslide and the smaller slip used in the back analysis. We applied a third-order moving average (Turmel et al. 2018) to determine the average retrogression distance along both the eastern and western banks of the rivers. Thus, the resulting estimate varies with location along the river. The landslides along the western bank showed greater retrogression than those found along the eastern bank (Fig. 12). Note that the 1971 landslide occurred on the western bank.
From this graph, we can map zones having a landslide risk, which takes into account the variable-width areas related to the retrogression distances of past landslides (see Fig. 13). The results of this method are more realistic, as they take into account the actual characteristics of the geotechnical and hydrological environments. To design this map, we considered all ≥ 14° slopes having a height of ≥ 4 m, following the MTMDET (2017). The widths of the areas classified as at risk of landslide were variable; we observed a constraint zone of approximately 80 m wide for the western bank and about 40 m wide for the eastern bank (Fig. 12).

Local slope stability analysis of the SJV area
The database, the 3D model, and the zonation map created for the SJV area facilitate further stability analysis of slopes across the area. A stratigraphic section can be created for a local study site. A back analysis can then be conducted at an existing landslide site to estimate the soil geotechnical parameters at a true scale and identify the appropriate analytical method to use for calculating slope stability. Finally, the mapping of the constraint zone will identify the main zones at risk of landslide, and any project proposed for the SJV area must first take into account this map of landslide risk (Fig. 13). The site of planned infrastructure located in a constraint area should be investigated further through a local stability analysis. The database could then be used to generate a 3D model of the site, and this model could then be imported into modeling software to assess slope stability. The data acquired from the back analysis of past landslides helps assign realistic values to the geotechnical parameters used in the stability analysis.
As an example, we conducted a local stability analysis of the slope near the site of the planned lookout (Fig. 14). Its projected location is at the top of an embankment characterized by a landslide risk. We obtained the topography of the section from a 3D model of the site that also indicates that the slope is composed entirely of clay.
Our modeling of the lookout slope, using the Slide software, assumed a uniform load weight, and we applied the "Corps of Engineers #2" method, as it was the most appropriate for this analysis on the basis of the earlier back analysis. The geotechnical soil property values were also determined through the back analysis: a fully Fig. 11 Map of stress zones related to landslides 22D06-050-0805 and 22D06-050-0705, 1:5000 (Ville de Saguenay 2007). 1, Aux Vases River; 2, Petit-Bras River; orange: embankments having a height ≥ 5 m and a slope of ≥ 14°, composed of soils of an undetermined nature with or without erosion, likely to be affected by landslides of natural or anthropogenic origin; pale yellow: 60 m protection zone at the top and foot of the slope saturated clay soil, γ = 18.6 kN/m 3 , φ = 35°, and cohesion = 8 kPa. The uniform weight of the lookout was estimated at 34.01 kPa. The modeling (Fig. 14) produced a minimum safety factor of 1.31.

Discussion
LiDAR data have been very useful for mapping the topography, increasing the density of the data points collected from a study area, and enhancing the precision of 3D modeling. Slope geometry is a crucial aspect of stability analysis, as it is a controlling factor of a possible slide, and the precision of LiDAR-derived data facilitates the detailing of slope geometry.
Adapting the electrical resistivity reference ranges to the study area was more difficult than expected. This difficulty can be partially explained by the complex structure of the Saint-Jean-Vianney soils. The clay is stratified and has been affected by at least two landslides; one partially remolded the soil and created inclined blocks, and the other landslide locally reshaped almost the entire soil mass. Our limited time for fieldwork also complicated the task of adapting the electrical resistivity reference ranges. Ideally, we would have included a greater number of control points. The TEM approach allowed us to obtain data at a sufficient depth; however, the spatial precision of the results was limited at less than a meter-scale resolution. For a specific point on the landscape, TEM results are certainly less accurate than borehole data. Borehole data is accurate to the nearest centimeter, while TEM data comes from surveys whose spacing becomes larger with depth.
Here the order of magnitude of the spacing can vary from decimetric to metric. Nevertheless, TEM is highly useful for interpreting the stratigraphy between boreholes.
The 3D model, as a tool for geotechnical investigation, must fulfill three important tasks for stability analysis: (1) accurately reproduce slope geometry; (2) estimate the depth of the bedrock; and (3) differentiate the major soil deposits, namely clay, sand (granular), and bedrock. The 3D model can be improved by using more exact parameter estimates. We simplified our categorization of the regional sedimentary layers to match the degree of precision of the TEM results. The topography of the rock underlying the study area, however, must be reassessed. The absence of rock data among the control points and the electrical resistivity values of the other interpreted sedimentary materials-moving downward through the accumulated material-suggest that the electrical resistivity values for the rock have been overestimated. The scattered occurrence of rock "peaks" reaching toward the surface, as observed in Figure D, reveals that the spacing between the TEM survey points was too great in some areas to reveal continuity between those "peaks". The back analysis procedure using limit equilibrium methods show that the computed failure surface did not superimpose perfectly on the actual interpreted surface. Nevertheless, the shapes of both surfaces were quite similar and were in close proximity to each other. This difference in the calculated locations of the failure surface-our estimated parameter values at failure were a cohesion of 8 kPa and a friction angle of 35°-likely relates to the accuracy of the pre-failure slope geometry. This geometry had to be assessed using topographic estimates derived from aerial photos collected between 1971 and 2015. We also identified that the "Corps of Engineer #2" stability analysis method provided appropriate results and can be used for further stability assessments for the region. We used an analytical approach because the estimated parameters are quick to compute, and this approach can produce sufficiently reliable results. More accurate results can be obtained using numerical software, which can model the progressive failure mechanism and are better suited for retrogressive landslides in sensitive clays.
Our mapping of the constraint zones followed largely the procedure used by the Saguenay city. We identified all slopes ≥ 14° having height of ≥ 5 m. Using the 3D model, we determined the main material composing each slope; this was clay in most cases. Along the Aux Vases River, we determined the constraint zones through the statistical analysis of previous slides. For other locations where there was no major vector of erosion, such as a river, we applied a 30 m wide constraint zone. The produced zonation map can determine whether a planned project lies within a constraint area and whether the development proposal requires further investigation. The new parameters determined by back analysis could be used for doing a more sophisticated analysis for these zones.

Conclusion
We presented a methodology for developing a landslide hazard assessment across a territory. For this approach, we proposed three tools: (1) a 3D model produced using a geophysical (TEM) survey and an adapted electrical resistivity chart; (2) the estimation of true-scale geotechnical parameters using a back analysis of past landslides in the area; and (3) a zonation map of slopes at risk of landslide, which includes constraint areas at the foot and top of these slopes. The construction of an appropriate 3D model demonstrated that our methodology is effective. We acquired approximately 75 new stratigraphic logs within a 12-day period, without using a drilling campaign. Results could have been improved even further if we had obtained a greater number of control points and more data related to bedrock depth, allowing further adaptation of the electrical resistivity chart. Overall, the developed 3D model demonstrated its utility as a tool for stability analysis.
Back analysis provided estimates of the geotechnical parameters at the slope scale. A cohesion value of 8 kPa and a friction angle of 35° produced conclusive results with the stability modeling of a slope affected by a previous landslide. We found the U.S. "Corps of Engineers #2" analytical method to be the most reliable approach for our study area. Nonetheless, a numerical modeling method that considers the progressive failure mechanism would provide a more exact output and would offer a more reliable approach for slope stability analysis in sensitive clays.
The zonation map constitutes an initial investigation tool for development projects within the study area. The 3D stratigraphy model, which served as the basis of the zonation map, can be combined with the parameter values obtained through other methods described in this paper to undertake stability analysis on other slopes within the study area.