Shear strength testing of consolidated claystones: breakpoint detection of shear stress versus shear displacement curves, a statistical approach

The evaluation of shear stress versus shear displacement curves is in the main focus of geotechnical engineering. Such curves, depending on the rock assessed, consist of a quasi-linear section, followed by a “kick” representing the peak shear strength, and a residual part, mostly parallel to the abscissa. The aim of the present study is to facilitate the future automatic detection of these crucial characteristics to take a step towards replacing their visual/analogue determination via modern statistical tools. Breakpoint detection methods (Cross-Entropy, Change Point Model) were applied to curves obtained from laboratory shear tests describing the shearing along discontinuities of nine Mont Terri Opalinus Claystone samples. Smooth and moderately rough claystone surfaces were studied. Results indicated that the end of the rising section and the kick observed on the shear strength curves was effectively approximated with the Change Point Model framework. An additional practical advantage of applying statistical tools such as breakpoint detection to shear strength determination is that it ensures the comparability of the obtained results.


Introduction
Shear strength is a crucial parameter in assessing the behavior of soil and rock masses in geotechnical engineering. Shear tests of consolidated soils and rocks along discontinuities provide essential input for both (i) the design of new surface-and 1 Page 2 of 18 subsurface engineering facilities (Crisci et al. 2019), and (ii) the stability calculation of natural formations in consolidated soil or jointed rock environment (Hencher and Richards 2015;Singh and Basu 2018), e.g. consolidated claystones. Besides the properties of the intact soft and well-cemented rocks (Contreras et al. 2018;Zhang et al. 2018), the properties of the actually jointed rock and rock discontinuities are also very important (Casagrande et al. 2018;Renaud et al. 2019). One such property is the shear strength along discontinuities, which is a function of the surface roughness, dilatancy and weathering of the rock surface (Barton 2013;Ulusay and Karakul 2016).
Previous analyses mainly dealt with the analysis of rock-rock (Grasseli and Egger 2003;Casagrande et al. 2018;Kim and Jeon 2019;Morad et al. 2019) and concrete-rock interfaces (Ghazvinian et al. 2012;Renaud et al. 2019). Approaches focusing on analyzing surface roughness (Muralha et al. 2014;Özvan et al. 2014), using photogrammetry (Casagrande et al. 2018;Kim and Jeon 2019;Morad et al. 2019) and statistical analyses of surfaces (Cai et al. 2018;Renaud et al. 2019) are also common. It is a fact that samples of different origin characterized by the same joint roughness coefficient can produce very different strengths. This implies, that shear stress versus shear displacement curves rarely constitute results that can easily be approximated with well-defined shear strength envelope curves (Hencher 2015). This creates a challenge for their objective evaluation, due to the fact that, without defining the peak and residual strength it is not possible to draw the envelope curve.
The application of statistical tools to better understand shear strength measurement results has not yet to become common, despite the increasing number of studies dealing with the topic (Casagrande et al. 2018;Morad et al. 2019;Renaud et al. 2019). In general practice, the break points (peak and residual strength) of shear strength-shear displacement curves are evaluated manually in a time consuming and quite subjective way which has numerous limitations (Muralha et al. 2014). The statistical detection of abrupt changes (breakpoints) in shear stress versus shear displacement diagrams has not yet been automatized. Previous works focused mostly on the effect of surface roughness on shearing (Grasselli and Egger 2003;Barton 2016;Cai et al. 2018) rather than the statistical evaluation of the shearing curves. In this paper we suggest a statistical toolset not applied before in shear stress versus shear displacement detection, which allows saving significant amount of time and resources and ensures comparability.
The Opalinus Clay is an excellent example of lithology to test shear strength along discontinuities, since its lithology is relatively homogeneous (Crisci et al. 2019;Nussbaum and Bossart 2008), due to its very small grain size (Favero et al. 2016), thus overcoming the problem of textural inhomogeneities caused by the presence of large crystals or a very heterogeneous mineralogical composition alongside a micro-fabric. The discontinuity surfaces present a non-uniform roughness; based on the type of irregularities they feature, previous studies suggest that two types of discontinuity surfaces exist: smooth and moderately rough (Buocz et al. 2014(Buocz et al. , 2017)-as explained also in Sect. 2.2.These factors allow the testing and validation of how effectively and accurately statistical methods can be employed to describe the different shear strength / shear displacement curve characteristics. The present study aims to suggest a new approach to evaluate shear stress versus shear displacement curves. It is designed to facilitate the future automation of the determination of the linear section(s) in the output of direct shear tests along discontinuities using various breakpoint detection methods. Via breakpoint analysis, the characteristic points of the stress-deformation curves are determined, and using the same methodology the results will become comparable between the evaluation of shear strength of different geo-materials.

Materials and methods
The tasks of the present research fell into two main parts: the preparation and processing of the samples to acquire the data, which included the direct shear strength tests in the laboratory; and the statistical analysis conducted on the data preprocessed for breakpoint analyses (Fig. 1). The Materials and Methods and Results sections were organized along the lines shown in the flowchart ( Fig. 1) describing the various steps of the analyses.

Geological background
The studied rock was Opalinus Claystone, which is exposed in the Jura Mountains of the Swiss Alps (Fig. 2). The Aalenian (Middle Jurassic) Opalinus Claystone forms a part of the Folded Jura system and is covered by Middle Jurassic limestone and Late Jurassic clay and limestone. The folding is related to the Alpine orogeny (Upper Miocene-Lower Pleistocene). The Opalinus Claystone is an overconsolidated clay (Favero et al. 2016), which has an apparent thickness of 160 m and, nowadays, an actual thickness of 90 m (Nussbaum and Bossart 2008). The nine samples studied are specifically from Mont Terri Underground Rock Laboratory and considered as shale due diagenetic processes that caused partial bonding of the clay. It is composed Fig. 2 The Opalinus Claystone beds and the location of Mont Terri Rock Laboratory in Jura Mountains (Switzerland) (taken from Nussbaum and Bossart 2008) of clay (66-67%) with a low amount of quartz (9-10%) (Favero et al. 2018). The obtained samples contained natural open discontinuities with no fill.

Laboratory tests: derivation of shear strength data
The Opalinus Claystone specimens were embedded in mortar for fixing in the shear apparatus. For all specimens, the bottom sample parts were digitalized together with their shear holder box using a 3D photogrammetrical method which allowed for obtaining a 3D point cloud of the samples ( Fig. 1: Step I; Buocz et al. 2017). The procedure consisted of three phases: 3D image production, referencing and data processing. The two images taken were matched together from the samples along matching points and the 3D image of ~ 30,000-60,000 points was created; for a detailed description on the calculations described above; see Buocz et al. 2018. The referencing was carried out with reference points on the sample, so that the software could calculate the real distances and orientation in a global coordinate system. The points of the sample surface were extracted from the point cloud. A linear regression plane was fitted to the point cloud of the sample surface ( Fig. 1: Step II; Fig. 3) with least squares fitting. This determined the shear plane (the linear regression plane of The determination of the surface roughness ( Fig. 1: Step III) was calculated as follows: i. Distances were calculated from the surface points to the regression plane. For minimizing failure due to the uneven distribution of the points in the point cloud, one single distance value was calculated in a 1 × 1 mm grid on the regression plane. ii. The surface roughness classes were calculated as frequency diagrams of the distance values between the linear regression plane and the surface points obtained for each 1 mm 2 cell. The frequency diagrams have a resolution of 0.1 mm ( Fig. 4).  Based on the analysis of the distances between the plane and the surface points, two categories of roughness were defined smooth, and moderately rough discontinuity surfaces (Buocz et al. 2017). The samples with smooth surfaces are the ones where the vertical amplitude of surface irregularities was between ± 0.5 mm, while samples with moderately rough surface were defined for larger values of this amplitude, in the range of ± 2 mm. In the smooth surfaces 50% of the cumulative frequency values are less than 0.25 mm distance from the regression plane and more than 90% are less than 0.5 mm (Fig. 4). The cumulative frequency values reach the median between 0.25 and 0.55 mm and over 85% of the distances from the regression plane are less than 1 mm (Fig. 4).
The area of the samples was measured digitally and later on used to calculate the area of the sheared surface ( Fig. 1: Step IV). The dimension of the samples was 5 × 7 cm. The vertical amplitude of the discontinuities varied between ± 2 mm. Next, direct shear strength tests were performed on nine pairs of rock samples with natural discontinuities, under a 1 MPa constant normal load ( Fig. 1: Step V) following the guidelines of the International Society for Rock Mechanics and Rock Engineering (ISRM; Muralha et al. 2014). A 1 MPa constant normal load was selected in order to detect even minor changes during the shearing tests. Three vertical and three horizontal displacements (two perpendicular and one parallel to shear directions) were measured using linear variable differential transformers adjusted to the upper shear box. An additional device was measuring the displacement of the bottom sample holder box moving in the direction of shear with a constant velocity of 0.8 mm min −1 (Fig. 1: Step V; Fig. 5).
The values of the shear load were recorded every 0.4 s throughout the course of the direct shear test. Calculation of shear stress and normal stress was performed based on the shear and normal load taking into account the changing surface area ( Fig. 1: Step VI). The obtained shear stress-shear displacement diagrams and the surface roughness classification served as the input of further analyses.

Break point detection: statistical methodology
Breakpoint detection aims to find abrupt changes and inhomogeneities in data streams of various fields of science. Most frequently, abrupt changes in the mean, variance and trend of the data are investigated. In earth sciences, the application of breakpoint detection techniques is most typical in climatology (e.g. He et al. 2008, Ruggieri et al. 2009, Lu et al. 2010, Topál et al. 2016Izsák and Szentimrey 2020), however, examples can be found in e.g. petrophysical research (Topál et al. 2016), geodesy (Gazeaux et al. 2015), geophysics (Topál et. al. 2017) or remote sensing of hydrological data (Bock et al. 2020). Despite the growing number of applications and the increasing number of scientific fields taking up this approach, there has not yet been any example in geotechnical-engineering and rock mechanics, specifically in aiding the evaluation of shear stress curves. Most of the studies focus on the question of the exact determination of the peak and residual shear strength by developing different models and formulae (Barton 2013(Barton , 2016Grasselli 2006; 1 Page 8 of 18 Liu et al. 2015;Casagrande et al. 2018;Shang et al. 2018) rather than detecting the peaks of the shear strength curves.
Since there is no "one-method-fits-all" to detect breakpoints and this is the first geotechnical-engineering application such methods were chosen which (i) are applicable to the statistical characteristics of the data at hand (ii) have been applied in various fields of science and (iii) have been tested on artificial data (Topál et al. 2016).
One of the methods used, is a multiple breakpoint detection approach, which is a model based on iterative stochastic optimization procedure relying on a modified Cross-entropy method (CEM; Rubinstein and Kroese 2013) employing a modified Bayesian Information Criterion to estimate the number and location of breakpoints in the mean and in the variance (for details see Priyadarshana and Sofronov 2015). With every iteration of the investigation in R statistical environment (R Core Team 2018), the parameters of the sampling distribution are updated using an elite sample (rho), until a stop criterion is reached (eps) (for details see Priyadarshana and Sofronov (2015). Given to the continuous nature of the data, the CE.Normal. MeanVar() function was used with distyp() set to 2, while the maximum number of breakpoints (Nmax) was set to five, and rho = 0.25. All other options were left as default. The parametrization of the breakpoint detection methods was chosen based on the comparative review of Topál et al. (2016). It has been observed that after each run, the CEM produces a slightly different result due to its stochastic characteristics, thus, in every case the number of runs was 100, and of the distribution of results, the three most abundantly occurring breakpoints were taken as the final result.
The other method, the Change Point Model (CPM) framework, was developed to detect change points in variance and/or the mean by evaluating the distribution of the data before and after a given point in the data stream using a two-sample hypothesis test (Shirvani 2015). It incorporates numerous tests specifically tuned for different types of statistical distributions, e.g. normal, (possibly unknown) non-Gaussian. Therefore, a breakpoint occurring at a particular time ( = k ) is detected if the distribution of the data is different in the segments before {X 1 , … , X k } and after {X k+1 , … , X n } for a finite sequence of independent random variables {X 1 , … , X n } , where x t is a particular realization of X i X i at a given point in time t = i (Shirvani 2015). The specific settings for the tests using the processStream() function were as follows: average run length (ARL0) was set to 200, corresponding to = 0.95 , the startup sequence (S) was n/1.5, where n denotes the length of the data stream.
Both methods (CEM and CPM) are capable of handling a stream of observations and multiple breakpoints, and it should be noted that based on experience (Topál et al. 2016) the locations of the breakpoints indicated by the tests give the locations of the last datum occurring in the sequence before a particular change takes place. For further information on the rationale behind the specific settings the reader is referred to the study of Topál et al. (2016) and the original documentation of the used methods referred to in Sect. 2.5.

Data preprocessing, prior to breakpoint detection
As a statistical-preprocessing step of the data prior to breakpoint detection, the first differences of the shear stress and shear displacement curves were determined to avoid any bias in the detection caused by the possibility of linear drift in the data, since the methods employed for change point detection assume identical distribution between the change points. Therefore, if the gradient of the linear slope is quasiconstant, it is to be expected that a change point will be found at the end of that linear section and/or around the peak shear strength corresponding to a change in the slope.
In order to be able to decide whether tests from the CPM framework for normal, or non-normal distribution should be used, the distribution of the shearing data had to be tested. This was done by Shapiro-Wilk test for normality (Shapiro and Wilk 1965) and by plotting Cullen-Frey diagrams (also known as Pearson plot) with bootstrapping (Cullen and Frey 1999) performed 5000 times in the present study to obtain an accurate picture on the distribution of the data. Latter was used to obtain a summary of the properties of a distribution of the data at hand.

Software used
3D photogrammetry was performed in ShapeMetrix3D; for a detailed description see Gaich et al. (2014) and/or Pötsch et al. (2019. The sample surface was estimated using the mvtnorm (Genz and Bretz, 2009), rgl (Adler et al. 2020) and scat-terplot3D (Ligges and Maechler 2003) packages in R (R_Core_Team 2018). The roughness of the surfaces was determined using a PHP script of the authors, while the surfaces were digitally mapped with AutoCAD 2016. The breakpoint detection was performed in R statistical environment (R Core Team 2018) using the Table 1 Coefficients (A and B), intercept (C)-all significant at p < 0.001-and fit statistics of the regression planes calculated for the shear planes of the studied Mont Terri Opalinus Claystone samples The surface types of the samples are indicated in the rightmost column. The surfaces of the samples marked with an asterisk (*) were partly dissected resulting in a poorer fit

Results and discussion
In previous works, surface calculation is performed using another approach in which it is recorded by two cameras integrated into the measurement head, from two different angles how various white-light fringe patterns are projected onto the object surface (Grasselli 2006), while the presented procedure provides new documentation of the shear plane presenting comparable and reproducible results. After approximating the surfaces and the shear plane of the nine samples, with linear regression, significant (p < 0.001) models were obtained categorizing the surfaces of the samples: smooth or moderately rough (Table 1). These results are shown on frequency diagrams (Fig. 4).
As an initial step in breakpoint detection, non-normal distribution was significantly (p < 0.0001) confirmed by the Shapiro-Wilk test for both the smooth (Fig. 6a) and moderately rough samples (Fig. 6b). In the analysis, the 5,000 bootstrapped values did not scatter but rather grouped together in the Cullen-Frey diagrams (Fig. 6), rendering the present result even more robust. This method is frequently applied in other fields, especially in systems (e.g. Nouri et al. 2014) and health science (e.g. Pouillot and Delignette-Muller 2010). It should be noted, however, that bootstrapping was not applied to those either. Therefore, thanks to the characteristics of the shearing data, the Mann-Whitney test was chosen and applied within a Fig. 7 Breakpoints detected on the S1-1 a, S8-1 b, S8-4 c smooth surface and on the S2-1 d, S2-2 e and S3-3 f moderately rough surface Mont Terri claystones. The secondary y-axis indicating the counts, refers to the number of occasions the CEM pointed a breakpoint to the given shear displacement out of its hundred iterations (red squares). Only the three most abundantly occurring breakpoints are marked in the figure regarding the results obtained with the CEM (color figure online) 1 Page 12 of 18 CPM framework (CPM-MW) tuned to detect changes in the mean of non-normally distributed data (Ross et al. 2011). Parallel to this, the application of the most frequently used tests within the CPM framework-to detect changes in the mean of normally distributed data (Student and GLR)-were omitted (Barton and Choubey 1977).
With the breakpoint analyses having been run, the MW test found a single breakpoint mirroring the peak shear strength in the case of the samples with smooth (Fig. 7a-c) and moderately rough (Fig. 7d-f) degrees of surface roughness. The test however, indicated two change points in only one case, the moderately rough sample where the peak area was elongated (Fig. 7f). Nevertheless, these two breakpoints clearly indicated the beginning and the end of the peak section.
The other breakpoint test method chosen, CEM found peaks in the curves, while indicating a breakpoint at the end of the linear "plateau" section after the peak. In the case of shear strength curves with a clearly concave rising section (Fig. 7c, d), it found a breakpoint at the inflexion points.
The curves describing the change of rock behavior under stress consist of subsections with different characteristics (Muralha et al. 2014). In the course of the shear strength tests, the one or multi-axis shear stress curves obtained have both linear and non-linear sections as seen in the case of the studied nine Opalinus Clay samples (Fig. 7). When breaking, peaks are formed, before and after which the shape of the curve (Barton 2013;Barton and Choubey 1977) and the variability of the values differ greatly (Fig. 7f). Muralha (1990) described different material models that are linked to the shape of shear strength versus shear displacement curves, with various different curve geometries strongly associated with bilinear, elasto-plastic, trilinear and curvilinear models. The determination of these characteristic points (e.g. the starting and ending points of the linear section, the peak shear stress etc.) is a crucial task in the process of shear stress curve examination (Barton 1973;Barton and Choubey 1977).
The detection of changes in the characteristics of shear stress and shear displacement curves is of increasing importance in the process of so-called multi-step analyses (Muralha 2012;Muralha et al. 2014), in which multiple loading levels are determined with the use of one sample (Morad et al. 2019). The load levels have to be changed prior to the break in order to proceed with the analysis of the sample without breaking it. In this case, multiple local maxima and inflexion points can be determined in one curve (Fig. 7). In practice, this is mostly done manually, based on the expertise of the investigator (Usefzadeh et al. 2013;Singh and Basu 2016), affecting the reliability and precision of the results. The number of such studies is vastly increasing, as well as the need for comparable results. Since it is primarily analog procedure comparability is not ensured, however, with the presented methodology this could be improved.
The characteristics of the shear strength/displacement curves depend on the surface roughness (Vattai and Rozgonyi-Boissinot 2020) and the rock material itself (Muralha et al. 2014). It is an essential parameter, especially in radioactive waste disposal site design, where anisotropic excavation damage zones are found. Opalinus Clay samples, with smooth and moderately rough surfaces, were investigated.
The location of the peak shear stress was clearly detectable with both applied breakpoint detection methods (CPM and CEM), as well as the inflexion point of the shear strength/displacement curves (Fig. 7). Local breakpoints observed on stress-deformation curves (Vásárhelyi 1999) represent minor losses in strength (e.g. Fig. 7c first breakpoint detected by CEM), such as the opening up of micro-cracks before the occurrence of major failure (Costin 1983;Cai et al. 2004). These minor peaks forecast major deformations and precede ultimate failure, as it has been shown in uniaxial strength tests (Martin 1993;Ferrero et al. 2008), as well as in triaxial strength tests (Cai 2010). Peak detection along discontinuities in shear strength tests has another importance since it refers to changes in the surface, namely, the reduction of surface roughness (Renaud et al. 2019) and the smoothing of the sheared surface. The detected peaks can serve as an auxiliary parameter for the training of artificial neural networks in modelling shear strength e.g. Moayedi et al. 2019. In the case of the Mont Terri Claystone, low shear stresses are recorded. Therefore it is a more significant challenge to detect breakpoints compared to high shear stresses of gneiss, migmatites (Renaud et (Singh and Basu 2018) or even to mortar (Ghazvinian et al. 2012). This is linked to the fact that the shear surfaces of claystones are much smoother, therefore it is more difficult to observe surface changes related to shearing. In turn, change in surface roughness lead to a change in shear strength along the discontinuity, as reduced surface roughness causes reduced shear strength (Barton 2013(Barton , 2016 The breakpoint detection methodology suggested here may also be employed in other mechanical tests such as the detection of failure at uniaxial compressive strength and in multiaxial strength testing (stress-deformation curves). Additionally, it might be very useful in the detection of the modulus of elasticity, since it is able to delineate the linear section of the stress-deformation curve automatically.

Conclusions
The study presents a new direction in facilitating the challenging practice of evaluating shear strength curves along discontinuities on the example of nine pairs of Opalinus Claystone samples. Measurements were carried out, producing different shear strength curves for smooth and moderately rough surfaces, each with characteristically different peak shear strength. Two different types of breakpoint detection methods were used to automatically detect these peaks while ensuring more robust results and reproducibility. On the one hand, the Change Point Model's Mann-Whitney method (CPM-MW) was employed to approximate the end of the linear sections; it was found that the peak shear strength coincided precisely with the first differential of the shear stress and shear displacement curves, though this method proved to be less robust in finding the peak associated to the major shear displacements. On the other hand, it appears that the Cross-Entropy method was less successful in assisting in this procedure, so it would seem that further rigorous and systematic investigations are required. The determination of these peaks is one of the main challenges to overcome before the shear strength along discontinuities may be quantified and the procedure can be automated in the future. Compared to manual detection, the suggested approach is a faster and more accurate way of identification of peaks, and thus recording of shear strength; more importantly, it ensures comparability across various lithotypes and studies. If the suggested procedure can then be automated, it would (i) allow the rapid evaluation of similar data, and/or (ii) lay the foundation of the continuous detection of the pre-break state in the course of the multi-stage tests during laboratory testing.