1. Introduction
The lowlands of southern Quebec are covered by marine clays deposited about 9500 to 8500 y BP following deglaciation and the marine incursion into isostatically depressed areas [
1]. The Laflamme Sea covered the Saguenay-Lac-Saint-Jean (SLSJ) region, Quebec, and deposited thick deposits of sensitive marine clays. These sensitive clays have a similar depositional and geological history to late-glacial marine clays deposited in Scandinavia. The specific feature of marine postglacial clays is exceptional sensitivity (ratio of undrained shear strength to remolded shear strength) to disturbances such as overloading, leaching, and earthquakes, whereby all the criteria are met to behave as quick clay [
2]. The presence of sensitive clay deposits and their ability to rapidly disintegrate from solid-state to liquid-state makes the SLSJ region an unstable area prone to landslides. Over the past few decades, the SLSJ region has experienced several landslides linked to heavy rains combined with snowmelt, some having caused loss of human life and property [
3,
4]. Torrential rains in July 1996 produced more than 1000 landslides in the region [
5]. These events have underscored the necessity of a better understanding of the geotechnical behavior of the Laflamme clays.
Several studies have assessed the characteristics of these marine clays to identify the factors controlling this geotechnical phenomenon. In the SLSJ region, till covers 40% of the area, followed by about 30% glaciomarine clays, and nearly 15% coarse postglacial sediments, which leads to much uncertainty in the development of correlations between various geotechnical parameters along with the depth. Bouchard et al. [
3] investigated the characteristics of these sensitive clays and determined that the observed pre-consolidation values originated mainly through preloading and erosion. Locat et al. [
4] examined the relationship between mineralogy and the physicochemical behavior of the Laflamme clay and other marine clays in eastern Canada. They found that an increase in phyllosilicates and amorphous materials contributed to increasing the liquid limit (ω
L), plastic limit (ω
p), and plasticity index (I
p). The Laflamme clays exhibit a very high sensitivity value, often exceeding 1000 [
6] and are less plastic with a greater fraction of silt particles (and a smaller clay fraction) than Champlain Sea marine clays deposited in the St. Lawrence Lowlands [
7]. The mineral content of these marine clays in southern Quebec has a thermal conductivity (k
s) range from 2.0 to 3.7 W/mK, a value mainly dependent on the quartz fraction (0.0–0.4) [
8]. These previous studies assumed that physical and mechanical parameters defined the properties of soils, which can be assessed through a classification system and the estimation of these parameters from collected soil samples. Clay deposits generally exhibit significant variability in geotechnical properties even within a seemingly homogeneous stratum as a result of unpredictable geological processes and assessing conditions such as field and laboratory testing. However, determining the geotechnical parameters through laboratory or in situ tests is always costly, and, depending on the location or the geological formation of the collected samples, uncertainties associated with geotechnical parameters are common. The uncertainty that arises from this inherent variability of clay properties greatly affects the mass displacement phenomena such as landslides and mudflows.
The uncertainties associated with a geotechnical parameter relate to their statistical distributions and the associated parameters, e.g., mean, standard deviation, and extreme values (
Figure 1). Outliers are defined by Barnett and Lewis [
9] as “observations which seem to deviate significantly from the other observations of the population from which they come, these observations seem to be inconsistent with the rest of the data, in relation with a supposedly known model”. Outliers can be random (stochastic) when they come from the inherent variability of the parameter or determined (epistemic) when they are the result of execution errors. Indeed, in
Figure 1a, neither the value x
(1) nor x
(n) seems to correspond to an outlier. In contrast, in
Figure 1b, x
(n) is an upper outlier located at the level of the right tail of the distribution. Contaminants are data originating from another sample but that are found inside the representative distribution of the currently tested sample; therefore, they are difficult to detect. They generally have a low impact on the statistical analyses unless the contaminant values lie at the extremities of the data distribution.
Figure 1c shows two contaminants indicated by a black circle; the one located on right is the extreme superior while that on the left is in the middle of the sample. Nevertheless, x
(n), although it is extreme and contaminant, is not an outlier. Finally, at the level of
Figure 1d, the extreme value x
(n), corresponds to a contaminant which is also an outlier. An outlier can therefore be the manifestation of the presence of a contaminant.
This paper presents a methodology to assess the uncertainties of geotechnical parameters at a regional scale. Being able to distinguish the type of outlier allows highlighting the impact of outliers in determining the uncertainties of eastern Canadian clays on a regional scale. This distinction then permits assigning a probability distribution to all studied geotechnical parameters. Several outlier treatments exist. Subsequently, assigning a probability distribution function to each parameter makes it possible to simulate the characteristics of clays with a given level of error. Different methods of treating outliers make it possible to observe the incidence of outliers on the Kolmogorov–Smirnov adequacy test [
10]. This methodology extracts the geotechnical information of clays without requiring any laboratory or in situ testing. Here we develop several steps to determine the uncertainties of geotechnical parameters at a regional scale. First, we detect outliers in the dataset, which has a considerable impact on the distribution type and the uncertainties of a parameter. After treating the outliers, we assign appropriate distribution functions to each parameter. We then produce correlations between selected geotechnical parameters for specific soils recovered from the SLSJ region. Finally, we discuss the selected parameters and their correlations.
2. Geology of the Study Area
The Saguenay-Lac-Saint-Jean region lies over the Grenville Province of the Precambrian Canadian Shield. The marine clays deposited in the Laflamme Sea 8500 y BP consist of pure clay as well as alternating bands of clay and sand [
12]. The visual description of the soil samples reveals the heterogeneous nature of grain size. Other surficial deposits in the region include sand and gravel, glacial till, and organic deposits (
Figure 2).
Our sampling sites lie in the municipalities of Chambord, Chicoutimi, Chicoutimi-Nord, Hébertville-Station, Jonquière, La Baie, Laterrière, Shipshaw, and Canton-Tremblay (
Figure 3). Soil samples lacking details of collection locations are designated as “unknown localities” or “unknown.” The data were obtained from laboratory tests carried out at the University du Quebec à Chicoutimi and from geological surveys of the PACES project [
13]. We developed a database of 6646 test results, including physical and mechanical parameters.
Figure 4 depicts an overall stratigraphy of the superficial geological units in the Saguenay-Lac-Saint-Jean region, Quebec [
14].
3. Methodology
We developed a methodology to characterize the geotechnical parameters of the SLSJ region, i.e., their statistical distribution, their uncertainties, and their correlations. This methodology allows us to (1) treat all outliers and to perform adequacy tests to assign an appropriate distribution law to each considered geotechnical parameter and thus determine uncertainties and (2) determine correlations between the various geotechnical parameters; we specifically studied arbitrarily chosen correlations to observe the influence of outliers on correlations.
3.1. Statistical Methods–Treatment of Outliers
We treated the geotechnical parameters individually and then grouped them according to test and soil type, regardless of sampling depth. We applied two outlier treatments: box-and-whisker plots (BM) [
15] and the Doerffel method (D) [
16], both following similar execution procedures. The Doerffel method requires estimating the threshold limit for outlier values (g), which depends mainly on sample size (
Figure 5). Doerffel’s method provides a graph for determining the threshold of outlier data values (g) for two levels of significance (Si) of 5% and 1% [
16]. To perform the Doerffel test, the average (
) and standard deviation of the data (S) are calculated regardless of the largest amount of data. Then, the largest amount of data (
) is considered to be outside of the row if it is true in the following equation:
where “g” is the threshold limit for outlier values which can be calculated by graph shown in
Figure 5.
The outliers must be detected in parallel using the equation specific to each of the two methods (
Table 1) for the data of each physical and mechanical parameter before being deleted. The identified outliers are corrected to be reinserted into the data distribution from which they were first extracted. Outlier processing produced four datasets as well as the initial data for each parameter. These datasets were (1) data with outliers removed by box-and-whisker plots, data with outliers removed by the Doerffel method, data with outliers modified by box-and-whisker plots, and data with outliers modified by the Doerffel method. In addition to the initial data, a geotechnical parameter can therefore have a maximum of five different databases depending on the existence or not of outliers in the initial data.
3.2. Kolmogorov–Smirnov Goodness-of-Fit Test
Determining the statistical parameters of a variable depends on the distribution law used to calculate these statistical parameters. The adequacy test is a process that discerns the adequate probability law for the descriptive analysis of a variable. We selected the Kolmogorov–Smirnov for the adequacy tests because unlike most other methods, this test can be used with any distribution law.
This test consists of determining the maximum difference between the cumulative frequencies of the observed data and the theoretical data, Dmax, a critical value, Dα, is determined based on the confidence interval and the size of the sample observed. The critical values can be estimated from Kolmogorov–Smirnov goodness-of-fit test (1) for a sample size of n > 80 and risk 5% by the formula 1.36/√n and (2) for a sample size of n > 80 and risk 1% by the formula 1.63/√n.
3.3. Determination of Uncertainties Associated with Parameters
The statistical parameters of mean and standard deviation can serve to determine the distribution of a variable. A geotechnical parameter containing outliers can be described by five datasets, each then subjected to an adequacy test. We, therefore, required a selection process (
Figure 6) following the adequacy tests to select the dataset for estimating the inherent uncertainties of this parameter. We computed the coefficient
Ck (
) for each dataset. The dataset with the smallest value of
Ck was used to determine the uncertainties of the geotechnical parameter relying on the probability law providing the best fit to this dataset.
The methodology of assessment of uncertainties of geotechnical parameters is detailed with liquidity index as an example.
- (a)
The initial data collected for the Liquidity Index (IL) was evaluated for outliers using D method and BM method. Liquidity index (IL) contained 5.8% outliers according to the BM method and 2.6% according to the D method.
- (b)
The outliers were modified in the next step according to the BM methods and D histograms were assessed. This process resulted in no visible outlier.
- (c)
If a few outliers are visible even after removal of the outliers using both methods, it was assumed that the appearance of new outliers is due to the inherent variability of the parameter.
- (d)
Use of underlying databases of outlier suppression is considered unwise when performing the adequacy test as this may alter an essential part of the information contained in the data.
- (e)
The other databases (BD whose outliers have been corrected) provided a conclusive result to the adequacy test. In the next step, the BD which offers the best adequate probability law for estimating the uncertainties linked to the liquidity index was determined.
- (f)
Further, the coefficient Ck () was calculated for each of the databases submitted to the test, and the uncertainties of the BD having the smallest Ck was representative of the physical parameter, IL.
3.4. Linear Correlations
We selected the geotechnical parameters to be analyzed based on correlations identified in the literature. We aimed to determine the relational strength between pairs of geotechnical parameters at both the regional and local scales. We analyzed the correlations according to whether the outliers were treated to identify the impact of outliers on the correlations. We assessed the correlation coefficient, the shape and equation of the correlation curve, and the data distribution of the two variables were projected onto scatter plots. SPSS V.25 software [
17] was used for these analyses.
4. Results and Discussion
The treatment of outliers does not have the same effect on the results of the adequacy test. The processing of outliers results in the parameters plastic limit, cone resistance, and sensitivity having no influence because it is the initial data that provide the best results. The assumption made with respect to the outliers detected in these cases is that they arise from the inherent variability of the geotechnical parameters in question. This assumption is valid as the sensitivity of marine clays can be very high. To verify the inherent variability of the other parameters, we referred to the published literature [
1,
2,
3,
18,
19,
20]. Our results for plastic limit and cone resistance matched published uncertainties. Moreover, we observed very high coefficients of variation (COV, >55%) and higher than for these three geotechnical parameters, thereby confirming the hypothesis that the outliers arise from the inherent variability of these geotechnical parameters.
On the other hand, we noted that the transformation of a geotechnical parameter can create or conceal outliers. The porosity data are proportional to those of the void ratio. Ideally, the distribution of void ratio should not have had outliers analogous to that of porosity. The concealing or creation of outliers risks distorting a descriptive analysis. The creation of outliers could, according to the operator’s judgment, require treatment of false outliers (substitution or deletion), and this judgment could bias the expected results. Concealed outliers would “contaminate” the distribution and lead to reliable results that are actually erroneous. Therefore, a perfect judgment for the transformation of a geotechnical parameter is required to identify the true source of any outlier as well as the authenticity of the transformation of the initial variable. In our case, we observed the best result in the adequacy test for void ratio data in which outliers were suppressed. For the geotechnical parameters studied,
Table 2 presents a summary of the treatment of VA determined using two methods, that of the box-and-whisker plot and that of Doerffel.
4.1. Uncertainty of Soil Physical Parameters
Uncertainties of soil physical parameters are certainly due to the inherent variability of soil that arises from unpredictable geological processes. Collectively, the origin of the uncertainty of geotechnical parameters initiates from the process of variability in depositional conditions and in the unpredictability of assessing conditions such as field and laboratory testing.
Table 3 summarizes the uncertainties obtained after the treatment of outliers and the adequacy test for all physical and mechanical parameters.
4.1.1. Natural Water Content and Atterberg’s Limits
For the water content of regional clays, we relied on 659 laboratory tests performed at UQAC. The data followed a normal distribution (
Figure 7) with an average of 41.79% (SD, 7.95%; COV, 19%). We observed the best results via the BM method when the outliers (VA) were removed. This suggested that the VAs were of epistemic origin. The obtained uncertainties gave values between 21.81% and 62%, which are much lower than published values for water content. For example, samples from the Saguenay Fjord had water content values between 50% and 154% [
4]. According to Sobotka [
18], the high values of water content are due to “a flocculated structure linked to a high void ratio”; however, after laboratory tests on samples having a void ratio between 1 and 2.58, the water content did not exceed 75%. This relatively low value for water content relative to the large index values of the voids implies a significant amount of air trapped in the interstitial voids of the clay samples.
We detected no outliers among the liquid limit data; the inverse Gaussian law was most applicable to the distribution of the data (
Figure 7). The slightly right-skewed distribution indicated that for most samples, the threshold of water content to pass from a liquid state to plastic consistency was relatively small, this content averaging 37.28%. Our data showed a low variability (COV 3%) and uncertainties varied between 16.5% and 67.11%. These SLSJ liquid limit values were similar to other regional sites, e.g., Leroueil et al. [
1] for Laflamme Sea clays and Windisch and Yong [
19] for eastern Canada clays. On the other hand, larger values have been recorded for marine clays in Scandinavia (73% with a SD 25.37%, Larsson 1980 [
21]), Finland (66.28% with a SD 19.75%, [
22]), and Sweden (71.06% with a SD 28.13%, [
22]). D’Ignazio et al. [
22] recorded values of around 200% for these Swedish clays.
Our plastic limits analyses relied on 160 regional sites. Only the box-and-whisker plot method detected outliers. The adequacy tests indicated a deterministic origin for the outliers; the outliers were detected by the box-and-whisker plot method because, despite treating the outliers, the best results using the adequacy test were obtained using the original data and Gumbel curve (COV 14%). Our data are again right-skewed, and our obtained values (11–28.4%) were similar to those recorded from other marine clays in eastern Canada and in Scandinavia, e.g., liquid limit of 26.35% recorded by Windisch and Yong [
19]. Clays from certain localities in eastern Canada, such as Kinburn, Quebec, produce plastic limit values greater than 30% [
20]. From the sample data available during analysis, we noted that above a water content of around 21%, a solid consistency would not be possible.
Our estimate of the plasticity index relied on 159 samples. We detected no outlier values. The Gamma law produced the best fit to the data, although variability was quite high (COV 61%; mean 16.9%, and SD 9.9%). Therefore, the quantity of clay minerals contained in the samples varies markedly with plasticity index values that do not exceed 44%. Our uncertainties obtained for the plasticity index produced lower values than those recorded for the Scandinavian clays [
19,
21]. Plasticity indices of 240%, for example, have been observed for some Danish clays [
23]. The clays in the SLSJ region appear less plastic than those in Scandinavia.
Despite the treatment of outlier values, high values for the liquidity index, e.g., 6.61, were preserved and showed much variability (COV 78%). Gumbel’s law provides the best fit. For the studied samples, more than 67% were in a liquid state and about 32% in a plastic state; a single sample presented a solid state (liquidity index of −0.69). Our obtained liquidity index values were higher than those found in the literature for clays from eastern Canada, e.g., 1.2 and 3.6 for clay samples from the Saguenay Fjord [
1], and for Scandinavian clays, e.g., around 1.2 [
19]. Therefore, samples from the SLSJ region are more liquid than Scandinavian clays.
4.1.2. Void Ratio and Porosity
We applied oedometric consolidation test (25 samples) to determine the void ratio. Tests were not carried out simultaneously on the same samples. Porosity is calculated by the relation: n = e
0/(1 + e
0), where n denotes the porosity and e
0 the data of the void ratio. For the samples measured by the oedometric consolidation test, the Gamma law produced the best fit to the obtained void ratio values (
Figure 8). Erlang’s law fits best when the VAs is corrected (oedometric consolidation test) via the BM method. The distribution was right-skewed, with an average of around 1. This value varied slightly depending on the COV (26%); nonetheless, it indicated that the void volume was often greater than that of the grains, hence the possibility of poor cohesion between grains. The Doerffel method did not detect any outliers for the porosity data measured using the oedometric consolidation test.
4.2. Uncertainty of Soil Mechanical Parameters
4.2.1. CPT Test
We relied on 1128 samples for cone tip resistance (Qc) and 1130 samples for cone lateral friction resistance (f
s) from cone penetration tests (CPT). We did not apply the Doerffel outliers’ treatment [
16] because this method did not match with more than 1000 samples. Concerning both geotechnical parameters, the cancellation of outliers detected by the BM method did not eliminate all outliers; we, therefore, presume that the origin of outliers is stochastic. The BM approach benefitted the fit test of the Qc data; Gamma’s probability law provided the best fit for Qc, and the log-normal probability law was the best fit for the f
s data (
Figure 9).
Cone tip resistance and cone lateral friction evaluate load-bearing capacity, and both increase with ground capacity and with depth. We recorded values between 60 kPa and 3250 kPa for Qc (average 1271.96 kPa, COV 56%) and values from 8 kPa to 87.5 kPa for fs (average 37.93 kPa, COV 66%). Nader [
23], in a study of Leda marine clays, found values between 100 kPa and 8000 kPa for Qc and 100 kPa and 650 kPa for f
s.
4.2.2. Preconsolidation Pressure
The preconsolidation pressure is the largest load that a sample has borne during its geological history. We relied on 43 samples for determining preconsolidation pressures in the SLSJ region. Doerffel’s method identified one outlier, and Erlang provided the best fitting probability law. The regional preconsolidation pressures showed a low variability and a mean value of 227 kPa (COV 10%;
Figure 10). These uncertainty values were slightly less than those recorded for eastern Canadian marine clays (average 212 kPa, COV 56%; [
19]) and Saint-Jean-Vianney (values between 200 and 1000 kPa; [
3]). Our uncertainty values for the SLSJ region were significantly higher than those of Finnish clays (average 96.03 kPa, COV 51%; [
21]) and Swedish clays (average 79.62 kPa, COV 72.6%; [
21]).
The overconsolidation ratio (OCR) indicates the degree of overconsolidation. OCR is defined as the ratio of maximum effective stress to which the soil has been subjected in its stress history to the present effective stress. We used 29 SLSJ samples for determining OCR. Only one outlier was detected by BM. Correcting this outlier produced a better result than its removal. OCR varied significantly, and Pearson V best fit the data (
Figure 10). We determine that about 79% of the samples are overconsolidated, agreeing with the findings of Bouchard et al. [
3]. Leda clays studied by Nadar [
24] were overconsolidated (average OCR 3.5) as were clays measured in Finland and Sweden [
21].
4.2.3. Sensitivity
Sensitivity tests provide the ratio of the shear strength of an intact soil (S
u) to a remolded soil (S
ur), and sensitivity (S
t) is the ratio S
u/S
ur to quantify the ability of a soil to flow. We used the field vane test (S
uv, S
urv, S
tv) and the Swedish cone penetration test (S
uc, S
urc, S
tc). We used 49 samples for the Swedish cone penetration method, 352 samples to determine the shear strength of intact samples via the field vane test, 61 samples for the shear strength of the remolded samples, and 403 samples for sensitivity tests. Regardless of the applied test, Sur and St had COV exceeding 100% or even 500% for the St parameter (
Figure 11).
Sensitivity values measured with the field vane and the Swedish cone indicated that most samples studied were highly sensitive and 46% of the samples were extremely sensitive (exceeded 1500). Suv values varied moderately (COV 42%), whereas those of Suc varied more (COV 74%) with respective averages of 99.27 kPa and 72.16 kPa. Uncertainties from the same parameters for Scandinavian clays are not as high as those in eastern Canada, which confirms the uncertainties of this study. Indeed, Scandinavian clays owe their resistance to cementation. Regardless of region, these sensitivity values are very high for all these regions and testify to the marked instability of clay deposits.
4.2.4. Unconfined Compressive Strength Test
Unconfined compressive strength (UCS) reflects the strength mobilized by the sample because of a vertical load. We relied on 38 samples for these determining UCS data. Only the box-and-whisker plot method found three outliers. A Weibull distribution provided the best fit for the UCS data (
Figure 12) with an average 43.2 kPa and COV 69%. More than 68% of the samples were from soils showing high consistency. The uncertainties obtained for UCS were similar to those of shear strength presented in the previous section for the SLSJ region, even though they are not the same samples.
4.3. Correlations between Geotechnical Parameters
For these correlations, we used 159 marine clay samples, 177 clay and silt samples, and 43 silt samples obtained for the region. According to Windisch and Yong [
18], Laflamme Sea clays generally have low plasticity, and all eastern Canadian clays have average plasticity. In this study, only linear correlations were considered for purposes of interpretation according to the Casagrande chart, and also for the comparison with other results, the nonlinear correlations of I
P according to ω
L were not considered in the literature. The R
2 value of more than 0.9 indicates a strong correlation between the two physical parameters and is found for all the other types of soil studied. In this study, clay and silt, clay, and silt samples varied from low to medium in their plasticity. Linear regressions indicated that physical properties varied widely among samples and that some samples were organic (
Figure 13). The linear relationship between liquid limit and plasticity index can also be seen in
Figure 13. The same observation was made for Leda clays [
22]. Likewise, the minimal difference between the clays of Finland and Sweden is that they are described as being very plastic [
20]. For clay and silt samples, we observed a greater dispersion of data; however, the correlations remained robust.
D’Ignazio et al. [
20] relied on this graphical representation to determine the number of samples for which water content was above the liquid limit. A similar trend was observed for the majority of Swedish and Finnish clays’ water levels and the liquid limits of these Scandinavian clays were double those of the Laflamme Sea clays (
Figure 14a). The mean value and the variabilities of ω and ω
L parameters for silt and clay material for different regions of SLSJ are presented in
Figure 14. The linear relationship ω-ω
L for clay samples exists regionally for the communities of Chicoutimi, Chicoutimi North, Jonquière (
Figure 14b), and unknown localities (
Figure 14c). This correlation observed for silt samples at Jonquière, Chicoutimi, and La Baie (
Figure 14d,e). We can observe from these figures that the regression lines obtained for clay samples had steeper slopes than those of the silt samples.
L’Heureux [
25] indicated that large retrogressive failures may occur in eastern Canada for I
L larger than 1.2 and/or S
ur smaller than 1 kPa. Lebuis et al. [
26] indicated that for soil with S
t greater than 25 and S
ur smaller than 1 kPa, a retrogressive flow slide is a possibility. In light of this, a correlation between liquidity index and sensitivity was developed and presented in
Figure 15. R
2 value of less than 0.78 for different locations indicates a weak correlation between the two geotechnical parameters and is similar for all the other types of soil studied. Similarly, for other studied locations, the relationship between the physical and mechanical was found to be nonlinear, which is not considered by this study.
5. Summary and Conclusions
Detecting and processing outliers are most effective when their flexibility is adapted to the inherent variability of the sample studied. The robustness of the statistical methods used in the study is significantly affected by the outliers. The Doerffel method qualifies to treat outliers of inherent variability origin. It is important to pay attention to the transformation of variables when manipulating to avoid the appearance or concealment of outliers; the operator’s good judgment is, therefore, essential in decision-making. The heterogeneity of the samples collected may have been a source of inconsistency between our results and those obtained from studies of eastern Canada clays, but similarities have been found between studied samples and samples obtained from other studies. The derived COV values for the Saguenay-Lac-Saint-Jean region are significantly high in comparison to the literature with cone penetration test data being only the exception. Comparison between our results with and existing information in the literature review raised slight differences between eastern Canada and clays elsewhere, especially for parameters such as water content, plasticity and liquidity indices, sensitivity, and pre-consolidation constraint. The SLSJ clays have less plastic as compared to the other Canadian clays as well as Scandinavian clays. The natural water content and liquid limit of Scandinavian clays are almost two times more than that of SLSJ clays. The evaluated uncertainties of geotechnical parameters help geotechnical engineers to assess the risk and stability of geotechnical structures, particularly in relation to mass displacement phenomena such as landslides and mudflows in the SLSJ region.