Comparison between methods for determining the effective vertical yield stress of intermediate fine-grained soils

Effective vertical yield stress (σ′xy) is essential in accurately describing fine-grained soils’ mechanical properties and their behaviour under load over time. It helps assess settlements and stress history. In most constitutive models, this parameter indicates changes in the soil behaviour due to the development of recoverable and irrecoverable strains resulting from loading. The results of oedometric compression tests for 25 soil samples with a wide range of plasticity parameters were used for the investigation. The intermediate fine-grained soils comprised different proportions of clayey, silty and sandy fractions. An in-depth, two-staged statistical analysis was carried out to compare twelve methods of determining effective vertical yield stress, namely: Casagrande (CM), Pacheco Silva (PSM), Butterfield (BTM), Oikawa (OIM), Onitsuka (ONM), Boone (BM), Becker (WM), Morin (WPUVSM) Wang & Frost (DSEM), Tavenas (SEM), Senol (SLSEM), and Janbu (JM). It aimed to check the association of these methods and the consistency of the obtained results. Based on the difference analysis, the methods originated in the work approach (i.e. WM, WPUVSM, DSEM) and CM gave comparable σ′xy values. The methods utilised bi-logarithmic plots (i.e. BTM, OIM, ONM) received slightly greater or lesser σ′xy values than BM and JM. The remaining methods were characterised by medium to the high variability and were sensitive to even the slightest disturbances resulting from the procedure of determining σ′xy.

In the field of geotechnical and geological engineering, the behaviour of soil is commonly assessed by means of laboratory investigation combined with constitutive modelling.The predictions for the soil response under any arbitrary stress path require an appropriate model's mathematical formulation.The elasto-plastic theory assumes that the material response is partly reversible and partly irreversible in response to applied forces.This implies the possibility of strains decomposing into recoverable elastic strains (ε e ) and irrecoverable plastic strains (ε p ).The transition between the elastic and plastic region is called yielding 1 , whereby this term should be regarded as a gradual transition from the elastic region to the plastic region rather than yielding occurring as a single point between the regions 2 .
The inviscid (rate-independent) elastic and time-dependent plastic behaviour can be described using elasticviscoplastic models, which belong to the family of general stress-strain-time models 3 .Sekiguchi 4 categorized the elastic-viscoplastic models into two main classes.The first one constitutes models based on the concept of overstress and corresponding Perzyna's overstress theory.The second covers elastic-viscoplastic models based on the concept of a nonstationary flow surface 5 .By using the concept of overstress, the decomposition of strains can be extended, taking into account the term of time, as the elastic strains are time-independent, whereas the inelastic strains are time-dependent.However, in such an understanding, the mathematical formulation of the model is based on strain rates so that the total strain rate is additively composed of the elastic and viscoplastic strain rates.It is generally accepted that the elastic strain rate is assumed to obey the generalized Hooke's law.In contrast, the viscoplastic strain-rate is expected to follow the non-associated flow rule.The main difference between the overstress theory and general elastoplasticity is the definition of inelastic strains.In the overstress models, inelastic strains are not related to the stress history but to the current stress state only, while in elastoplasticity, inelastic strains are related to the stress rate.
Despite these theoretical frameworks, from the practical point of view the transition from elastic to plastic soil behaviour is usually interpreted as a clear change in the stress-strain response.By convention the effective vertical yield stress (σ′ xy ) also known as 'pre-consolidation stress (σ′ p ) is used to characterize the compressibility

OPEN
Faculty of Civil Engineering, Cracow University of Technology, Cracow, Poland.* email: bartlomiej.olek@pk.edu.pl of soil.Casagrande 6 defined preconsolidation pressure as the maximum past vertical effective stress applied to the soil and correlated it to the change in a slope of void ratio versus logarithm vertical effective stress curve 7 .The term 'pre-consolidation stress' may be inappropriate because it shows that only one mechanical factor causes the change in soil behaviour during loading.The visible change in soil response may be a result of many factors, including not only the past depositional and stress history but also changes in water content, chemical-physical reactions (weathering, cementation, recrystallization of minerals, ion exchange, modification of the adsorptive water layer or intermolecular attraction forces of clay minerals), cold welding, bonding due to long-term secondary compression, ageing, and other diagenetic factors [8][9][10] .Moreover, some observations indicated that soil might show a so-called preconsolidation pressure much greater than the existing effective stress without evident erosion in its geological history 11 , or the pressure may decrease with increased soil depth 12 .Therefore, following Gao 13 , the term 'vertical effective yield stress (σ′ xy ) will be used in this work as more fundamentally correct.The vertical effective yield stress points to the stress level corresponding to the onset of an irreversible mechanical behaviour of soil, and hence is useful for analysing and predicting settlement behavior, overconsolidation ratio (OCR), stress history, and short-term stability problems in fine-grained soils [14][15][16][17][18] .When σ′ xy is greater than the current effective stress, the stress-strain relation exhibits softening behaviour 19 , and the soil is in an over-consolidated state.In turn, strain hardening occurs when σ′ xy is equal to the current effective stress, and the soil in this case is in a normally-consolidated state.
The present work deals with the issue of determining the vertical effective yield stress in the intermediate soils consisting of different proportions of clay, silt and sand fractions using available methods.An in-depth, twostaged statistical analysis was carried out to compare twelve methods.Its purpose was to determine the similarity of these methods and the consistency of the obtained results.The main goal of the research was to investigate whether the values of σ′ xy determined by different methods are repeatable (or similar) and evaluate the rank of discrepancies between the obtained values.The intention was to indicate practical approaches that can be used interchangeably and are the least ambiguous.In this respect, the soils collected from the selected area located in Poland were studied through incremental oedometer tests.

Overview of existing methods for determining σ′ xy
Over the past decades, various empirical methods for determining the value of the vertical effective yield stress (σ′ xy ) have been developed.The most popular methods include the graphical interpretation of the oedometric compression curve described by the relationship between void ratio and effective vertical stress and their different modes.This paper focuses on twelve different methods, namely: Casagrande (CM), Pacheco Silva (PSM), Butterfield (BTM), Oikawa (OIM), Onitsuka (ONM), Boone (BM), Becker (WM), Morin (WPUVSM) Wang & Frost (DSEM), Tavenas (SEM), Senol (SLSEM), and Janbu (JM).Based on previous observations, the conclusion has been made that all these methods might give different results when considering a particular compression curve.In general, this is related to the subjective assessment of the inflection point of the soil compression curve or the influence of the scale of the graph.The quality of the collected and tested soil samples is also essential in correctly determining the preconsolidation pressure.To minimize subjectivity in determining the compression parameters, the authors used the open-source software pySigmaP 20 .

Semi-logarithmic approach
Arthur Casagrande 21 was probably the first researcher who used the semi-logarithmic approach to estimate σ′ xy .Although many researchers judged his pioneering method as inaccurate 22 , it is still remains a standard for comparison to other methods 23 and recommended in widely used standards (i.e.ISO 17,892-5, ASTM D2435).To determine σ′ xy , the inventor used a semi-logarithmic plot of the void ratio e versus the effective vertical stress.The methodology for determining σ′ xy using the Casagrande method (CM) is illustrated in Fig. 1a.The graphical construction involves determining the minimum radius of curvature point on the compression curve and drawing several additional lines.5][26][27] ).It should be mentioned that semi-logarithmic approach has become the target of further improvements and refinements.For example Dawidowski and Koolen 28 developed computerized version of CM.Moreover, Peck et al. 29 , for sensitive soils, simplified the procedure, limiting themselves to determining the point of inflection on the curve and the point of intersection of the tangent line at this point with the so-called initial void ratio line.
In turn, Pacheco Silva 22 introduced the method called Pacheco Silva method (PSM), accounting for independence from the drawing scale 30 .The method utilised the initial void ratio line and graphical construction to determine σ′ xy , as shown in Fig. 1b.Grozic et al. 31 have demonstrated that the determined σ′ xy by this method depends on locating the tangent to the virgin consolidation line.
Further, Boone 7 noted that the CM does not always give accurate results, especially for soils that do not show a well-defined change in curvature in the graph.He also recognized that the derived value of the initial void ratio e 0 from the oedometer is not a favourable approach due to the sample's disturbance issue.On the basis of the above, the so-called Bonne method (BM) has been proposed.The methodology for graphically determining σ′ xy using BM is illustrated in Fig. 1f.The method requires establishing the in-situ void ratio e v0 corresponding to the in-situ effective vertical stress σ′ v0 (directly using load increment corresponding to σ' v0 or using the interpolation).The reader is referred to the original paper for a detailed empirical investigation and graphical procedure 7 .

Bi-logarithmic approach
Butterfield 32 has proposed an alternate approach for describing soil volume change with variations in effective vertical stress.It has been considered to plot ln(1 + e) versus ln mean effective stress (p′).The so-called bi-logarithmic method (BLM) used to characterise soil compressibility implies some advantages over the classic logarithmic plot in that it might be applied to large-strain deformation and clarifies the influence of the unloading stress ratio on reloading response and compression parameter values 33 .Moreover, the bi-logarithmic relationship might show improved linearity with experimental data for unloading and reloading sections of the compression curve.
In the case of determining σ′ xy , the linear portions at both ends of the original compressibility curve are extended.The point of their intersection (see point A) determines the sought value of σ′ xy .The bi-logarithmic nature of the compression behaviour has been reported for many different soils 9,18,[34][35][36][37][38][39][40][41] .The bi-logarithmic approach includes several methods for which the procedure for determining σ′ xy is the same, and only the form www.nature.com/scientificreports/ of mathematical relationship used differs.For instance, Onitsuka et al. 39 supported the use of the ln(1 + e) -log p′ plot, while Oikawa 42 , Jose et al. 23 and Sridharan et al. 9 preferred the log (1 + e) − log p′ plot.From the methodological aspect, the proposal of Onitsuka et al. 39 has strong theoretical background because it might be related to the work approach 43 .The methodology for determining σ′ xy using the methods of Butterfield (BTM), Oikawa (OIM), and Onitsuka et al. (ONM) methods are illustrated in Fig. 1c,d,e.

Work approach
Another interpretation approach has been proposed by Becker et al. 44 .The idea behind it was based on the definition of work (W) (total strain energy or energy per unit volume).Here this method is called as work method (WM).The authors defined incremental work (cumulative total strain energy) using the following formula: where: σ ′ i+1 + σ ′ i -effective stresses at the end of load increments i + 1 and I, ε i+1 , ε i -natural strain at the end of load increment i + 1 and given i-th.
Principles of WM are shown in Fig. 1g.The plot in an arithmetic scale in which the relationship between work per unit volume (W) and consolidation pressure (p) is expressed in bi-linear lines has been used to determine σ′ xy .The intersection of the so-called preyield (see AB line) and post-yield (see CD line) lines when the work incremental rate ΔW/Δp increases significantly has been interpreted as σ′ xy .
Morin 45 undermined the original way of computing the work per unit volume using cumulating incremental work values per unit volume.As an alternative, the work per unit volume of solids had been recommended (say work per unit volume of solids method (WPUVSM)).Hence, the two methods differ only in the mathematical formulation of work per unit volume for a given load increment, while the procedure for determining the σ′ xy is the same.
More recently, Wang and Frost 46 have improved the original WM.They considered the consolidation process irreversible, in which most of the work done during it is dissipated.Based on this foundation, the dissipated strain energy method (DSEM) was introduced to describe pre and post-yield behaviour based on plastic deformations.The schematic of DSEM in the specific strain energy (SSE) versus effective vertical stress space is shown in Fig. 1h.The determination of the σ′ xy does not differ from other methods based on the work approach except that it applies to the dissipated energy.

Other methods
So far, alternatively, to the semi-logarithmic, bi-logarithmic and work-based methods, some other methods of determining σ′ xy have been developed.Likely to the methods mentioned in the previous section, they are mainly based on observing strains exhibited by the soil material during oedometric tests.Certain exceptions are methods developed based on the time resistance approach 47 and the analysis of the variation of the tangent modulus (M = Δσ/Δε) as a function of stress.Unfortunately, the original description of the Janbu procedure (say Janbu method (JM)) is very vague.In some cases, the authors indicated the use of the relationship between M and consolidation pressure on a linear scale (see Fig. 1j), and in others, such as in the case of Jacobsen 48 , the relationship between deformation and consolidation pressure on a semi-logarithmic scale.In the first case, σ′ xy was obtained by locating a downward concavity on the graph.In the second, the experimental data was approximated by two straight lines, and their intersection pointed to the σ′ xy , similar to bi-logarithmic methods.In addition, two variants of the Janbu approach can be found in the literature (i.e 49,50 ).Other examples are methods designed for CRS studies 51,52 or based on shear wave velocity at small strains 53 .
Alternatively, researchers have developed methods based on the strain energy concept (SE = σ′ΔH/H).In this connection, Tavenas et al. 19 considered the intersection of two inclined lines on the plot of SE versus consolidation pressure in a linear scale as a σ′ xy (see Fig. 1k).Hereafter, in this paper, this method is denoted as the strain energy method (SEM).Moreover, Senol and Saglamer 54 studied the semi-logarithmic relationship between SE and consolidation pressure (see Fig. 1l).Thus, this method is denoted as the semi-logarithmic strain energy method (SLSEM).
To sum up, the past studies have all highlighted issues related to the practicalities of the graphical interpretation of consolidation test data.For instance, Schmertmann 27 , Brumund et al. 24 , and Holtz and Kovacs 26 pointed out that the semi-logarithmic approach applied to disturbed soils may be inaccurate due to the difficulty in determining the inflection point on the compression curve.Grozic et al. 31,43 , Clementino 30 and Boone 7 provided detailed comparisons between various methods.The main conclusion drawn from their considerations is that most of them are burdened with uncertainty and, in many cases, challenging and ambiguous to apply.To put it briefly, the technicalities of most of the methods are concerned with defining the tangent or "best-fit" lines and reading of yield point.The main error sources come from subjectivity and the plot type used for interpreting test results (i.e.scale effects).

Test materials
The soils used for the investigation were collected from the area that lies on the border of two polish geological units, i.e. the Carpathian Foredeep and the Outer Carpathians (see Fig. 2).On the basis of preliminary geological engineering analysis of the area, most of the sediments are Quaternary.In turn, the Miocene formations are the immediate bedrock of the Quaternary sediments.The samples were collected from the pseudo-loess cover.
(1) The core samples used for the laboratory investigation were taken from a depth of 1.65-7.85m.The cover is a typical example of the subsoil representing the area of interest.The accumulation of a substantial part of it is associated with the period of the northern Polish glaciations.The acquired soil material is mixtures of silt, clay and sand fractions in different proportions, mainly loams, silty loams and sandy loams with various interlayers, inserts, and smaller or larger clusters that differ in grain-size distribution, colour, water content and consistency state.The separated lithological-genetic types of soils, despite the general similarity in their formation, are more or less heterogeneous.When considering soil behaviour, these soils should be treated as "intermediate soils", which are neither clay-like nor sand-like.

Index properties and classification of soils
Comprehensive index tests were carried out to allege the soil identification, description and classification.These tests included determining natural water content, consistency limits, specific gravity, bulk density and unit weight.The gravimetric method with oven drying was used to measure the natural water content (w n ) as per ISO/TS 17,892-1 55 .The liquid limit (LL) and plastic limit (PL) were determined according to ISO/TS 17,892-6 56 and CEN ISO/TS 17,892-12 57 , respectively.Using the PL and LL values, the plasticity index (PI = LL-PL) (%) was calculated, which allowed for utilising Casagrande's chart for classification purposes.The particle-size distribution was identified prior to ISO-TS 17,892-4 58 , and the specific gravity tests were performed following the procedure described in ISO-TS 17,892-3 59 .It is worth noticing that organic matter in the tested soils was not found during the laboratory investigation.In the work presented herein, the soils were classified based on the Atterberg limits combined with criteria for delimiting clays, clayey soils and less plastic silty soils given by Moreno-Maroto 60 .

Consolidation tests
The compressibility characteristics of soils utilised in the study were determined by means of a conventional incremental loading oedometer test in accordance with the requirements of the ISO/TS 17,892-5 61 .The tests were performed on specimens placed in the consolidation ring with a size of 60 mm diameter and 20 mm thickness.The internal sides of the ring were lubricated with silicone grease to minimize side friction.Filter papers and porous stones covered each specimen's top and bottom faces.Testing soils were subjected to a loading-unloading path, and each load was kept for one day (24 h).The following loading scheme was executed: 12.5 kPa → 25 kPa → 50 kPa → 100 kPa → 200 kPa → 400 kPa → 800 kPa ← 12.5 kPa.The specimen deformation was expressed as relative compression (settlement) for each load increment, and the corresponding time was noted.Interpretation of the compression curves to determine σ′ xy was performed using an open-source pySigmaP software 20 .This made it possible to remove the interpreter's judgment when choosing inflection points or points of maximum curvature in the e − log(σ′ v ) space, necessary for methods belonging to the semi-logarithmic group.The pySigmaP, developed in Python 3, is designated to identify the yielding of soils by incremental loading oedometer test.This software enables to determination the σ′ xy of the soils using nine different methods presented above.The pySigmaP package is divided into six modules: data.py,casagrande.py,energy.py,bilog.py,pachecosilva.py, and boone.py.The data.py module contains the class to load and manage the data of the compressibility soil response.The remaining modules contain the classes and methods to determine σ′ xy .The pySigmaP provides a more objective and straightforward analysis to estimate σ′ xy .The undeniable advantage of the software used is its calculations are guided by analytical and numerical methods minimising bias by the scale dependency, analyst visual skills, engineering judgment, and experience.

Index properties and soil classification
Table 1 presents the index and selected physical properties of the tested soils.It may be seen that the values of natural water content ranged from 13 to 35%, LL from 27 to 69%, PI from 7 to 39% and the percent clay size fraction from 11 to 35%.Based on the values of PI and LL and the boundaries established by Moreno-Maroto, most of the tested soils are located between the C and M lines on the plasticity diagram (see Fig. 3).This indicates moderately or slightly clayey soils.The following 5 points are located above C-line, which shows the clay's limit.The remaining 4 points are below the M-line, which points to less-plastic silt-clay soils and mixtures of various fractions with the dominance of the silt fraction.As shown in Fig. 3, the classic Casagrande's A-line that separates clays from silts leads to utterly different classification results.Considering the visual inspection of the soils, their grain-size distribution and the values of inherent parameters, such as PL and LL, it should be stated that in most cases, we are not dealing with clays but with soils containing clayey fraction that only modifies engineering properties of that soils.

Compression of the soils
The semi-logarithmic plots of void ratio vs effective vertical stress and normalised void ratio vs effective vertical stress for 25 intact soil samples are shown in Fig. 4a,b.As can be seen, most of the e-log σ′ v relationships are initially curved concave downward.At the high values of the effective vertical stress, all exhibit well-defined straight lines or slightly pronounced concave upward course.The slopes of these lines are further utilised to derive the soil's compression index (C c ).
Inspection of the curves reveals some slight differences in the shape of their initial portions.Thus, the three patterns differ in the initial compression rate.The first case (type I) concerns mild changes in the parabolic shape of the curve with a not very clearly marked breakdown corresponding to a shift in soil behaviour (see soil OC13 in Fig. 5a).Therefore, the e-log σ′ v relationship is curved and concaves downwards throughout the stress ranges.This type is most common for the investigated soils (17 samples).It complies with the reported compression characteristics of compacted soils, soils with lower initial water content than their LL value 62 or less-plastic soils.The second type of curve (type II) is similar in shape to the first one, except for the initial segment, which is a distinct   rectilinear with a constant slope (see soil OC1 in Fig. 5a).This commonly may occur due to the significant soil disturbances created during sampling, transport and laboratory preparation/cutting of the specimen before the test.Compressibility curves classified as type II (i.e.OC1, OC15, OC16) in the plasticity diagram were located below the M-line, representing less-plastic silt-like soils.Finally, the initial curvature is flattened in the third case (type III), and the middle part of compression curve is convex upward (see soil OC25 in Fig. 5a).It should be noted that these differences are more visible when the vertical scale of the chart is reduced.Taking into account that the differences in the initial courses were minor and the middle and final sections of the compressibility curve are usually used for interpretation, the presented division into groups should be treated conventionally.
The compression behaviour of various fine-grained soils can be well represented by the two or three straight lines in the plot of ln(1 + e) against log σ′ v or log(1 + e) against log σ′ v 63-65 .These lines are used to represent pre-yield, transitional, and post-yield compression behaviour 66 .As indicated by Hong et al. 66 , the small-compressibility refers to the pre-yield regime (overconsolidated region, OC), where changes in the soil structure are relatively small up to the σ′ xy .The so-called "transitional regime" explained by transitional stress and gradual destructuration of soil structure can be observed in some soils.In turn, in the post-yield regime (normallyconsolidated region, NC), the soils experience the most remarkable changes in compressibility.Figure 5b presents oedometric inverse S-shaped compression data replotted in the graph of ln(1 + e) vs log σ′ v for the investigated soils.It is observed that the shape of all analysed compression curves can be imitated by the two straight lines.Thus, the soil's compression follows bi-linear behaviour.
Table 2 listed compression parameters for the tested soils.The values of the calculated C c for the stresses between 400 and 800 kPa were in the range of 0.057-0.252.In general, it can be stated that the higher the C c value, the greater the maximum straining of the sample, as depicted in Fig. 6.It is common practice to correlate C c with the Atterberg limits of the soil [67][68][69] the water content, the clay fraction CF (%) or the organic matter content.Due to the origin of the soil from different locations and the adopted utilitarian goals of the article, a detailed description of these relationships was omitted because the primary intention of the undertaken research is a statistical assessment of the methods for determining σ′ xy .In order to determine σ′ xy , the data from oedometric tests performed on intermediate soils were converted and replotted on the appropriate diagrams.The typical results of interpretation using the pySigmaP software for soil sample OC1 are shown in Fig. 7.The first step after loading the data was to plot the compression curves and determine their basic parameters, such as C c and C r .For the CM, a cubic spline function was used to determine the maximum point of curvature on the compression curve.As can be seen, this value was established as 202 kPa (see Fig. 7a).The use of the pySigmaP excludes the subjective interpretation of this point because the mathematical formulation for the curvature of the soil compression response as a function of the stress level defines it.Figure 7b shows the plot for the PSM.After determining the initial e 0 value, the elongation procedures of the linear part of the compression curve were performed.For the discussed soil sample OC1, the value of σ′ xy determined by the PSM was 63 kPa.The large scatter of values between these two semi-logarithmic methods   www.nature.com/scientificreports/ the relationship used, the same σ′ xy value of 138 kPa was obtained in these methods.Figure 7f shows the compressibility diagram for the BM.The value of σ′ xy , in this case, was 104 kPa and was smaller than those determined by the CM and bi-logarithmic methods.Figures 7g,h present the results of two methods based on the work approach for determining σ′ xy .Since both methods differ only in the mathematical formulation of the work per unit volume for each load increment, the determined values of σ′ xy for the analysed soil samples are the same.For the discussed soil sample, the σ′ xy was established as 147 kPa.This value was practically identical to the σ′ xy obtained with BM and close to the values determined by bi-logarithmic methods.Figure 7i depicts the diagram for the DSEM.Since this method uses dissipated strain energy that considers both elastic and plastic deformation, the obtained value of σ′ xy = 111 kPa was lower than those determined by WM and WPUVSM.The obtained value was also smaller when compared with other methods other than BM and PSM.
Figures 7j,k,l present the plots for determining σ′ xy by the JM, SEM and SLSEM.During the interpretation of the compression curve for the OC1 sample as well as for some other samples, problems with the application of these methods were encountered.In the case of JM, it was not always possible to locate the point where the curve bends downwards (i.e.method requires the judgement and experience of the interpreter).In addition, some relationships between M and σ′ v showed a practically increasing linear course.In the case of SEM and SLSEM, choosing a rectilinear segment and drawing tangents was problematic (see Fig. 7k,l).The results obtained with these two methods are ambiguous and difficult or even impossible to apply for many of the analysed compressibility curves.Moreover, with strict adherence to the guidelines provided by the authors of the SLSEM method, the determined σ′ xy values were significantly inflated compared to other methods.

Methods comparison study
In this subsection, a method comparison study is attempted.Due to the lack of a reference method admitted according to current practice as 'the best" and most reliable, the centre of attention was on assessing the degree of agreement between different methods and the precision of the methods.Method comparison is the statistical process by which error components are recognised and characterised.Thus, errors within acceptable limits may indicate the adequacy of the given method and point to its reliability.In this work two-step approach for comparison purposes was applied.At first, the scatter plots for analysing the comparative results have been constructed.The extreme outliers have been identified and, in some cases, excluded during the trimming of the data.The slope and intercept of the trimmed data using a regression procedure were also determined to consider the convergent validity of compared methods.In addition, the variance of the slope from the identity line that formed a 45° angle with the abscissa has been used to indicate a bias.The strength of the linear association between the two compared methods was evaluated by the Pearson product-moment correlation coefficient (r).In this way, the r-value helps establish convergent validity showing whether two methods are related to each other.Note that The calculated r-values were not used to measure the accuracy of methods or to assess the agreement between them.
In the case of compression laboratory data, neither provides an unequivocally correct determination of σ′ xy when comparing the two methods.Therefore, assessing the degree of agreement between these methods might be helpful.To evaluate the degree of agreement between the methods, the common practice is to study the mean difference between them and to establish limits of the agreement [71][72][73][74] .The limit of agreement is referred to the confidence interval that represents the range of values in which agreement between methods lies for approximately 95% of the sample.A difference plot and analysis of differences may be used to compare two "new" methods or one method against a reference one 75 .The second part of the analysis was based on establishing the statistical limits by using the mean and the standard deviation (SD) of the difference (indice of precision).In the difference plot, the Y-axis expresses the difference between the two paired values, and the X-axis represents the average.This allows for checking the assumption of normality of differences.

Association between the methods
A cross-correlation analysis was performed for the collected values of σ′ xy to assess whether the individual values determined by twelve different methods are statistically significantly related.Figure 8 presents the results of this analysis as a correlation matrix with scatter plots and computed r-values for the probability value p < 0.05.Due to the same procedure and nature of determining σ′ xy for three bi-logarithmic methods, a complete convergence of σ′ xy values was observed.The same relation was observed for the WM and WPUVSM.
The CM showed a low correlation with other methods of determining σ′ xy .The best link was found with the bi-logarithmic methods and with the SEM, for which the r = 0.693.In the same group of semi-logarithmic methods, σ′ xy using the CM correlated with those determined using the PSM with r = 0.688 and the BM with r = 0.649.The values of σ′ xy by the CM are associated with the SLSEM with r = 0.664 and the WM and WPUVSM with r = 0.639.The lowest correlations were obtained for the CM and DSEM, where r = 0.464 and the CM and the JM with r = 0.580.These results indicated that those correlations cannot be considered as significant.Many different factors can cause the low level of correlation of the CM with other methods.It is believed that most important is the deformation properties of the tested soil samples, for which the initial sections of the compression curve exhibited a relatively flat geometry and the accompanying scale effect.However, the subjectivity of finding the inflection point on the curve can be ruled out as the software determines it.It is also worth noticing that the CM gives higher values of σ′ xy for the vast majority of samples compared to other methods, which also affects the value of r.In the group of semi-logarithmic methods, the PSM and the BM showed a high correlation with r = 0.972.This may be due to the similar methodology for determining σ′ xy .The PSM is also well correlated with the SEM, for which r = 0.921, and the bi-logarithmic methods, for which the r was 0.923.
Satisfactory results were also obtained between the bi-logarithmic methods and the work methods.Concerning the WM and WPUVSM, the r values were the highest among all comparisons and amounted to 0.976.The r values for correlations between this method and the BM, the PSM and the WM and WPUVSM were respectively 0.910, 0.860 and 0.838.The values of σ′ xy obtained using the SEM correlated moderately with those from other methods, preferably with the PSM, for which the r was 0.921, and with the BM (r = 0.838).The association of the SEM and the bi-logarithmic methods, the WM together with WPUVSM and the DSEM were slightly weaker, with r amounting to 0.762, 0.755 and 0.738, respectively.The weakest correlations occurred between the SEM and CM with r = 0.693 and the JM with r = 0.507.The SLSEM correlated best with the SEM with r = 0.872 and PSM with r = 0.830 and was at the same moderate level of correlation.The BM achieved a value of r = 0.765, the WM and WPUVSM 0.755, the bi-logarithmic methods 0.693, and the DSEM 0.621.The SLSEM and the JM exhibited the weakest correlation with r = 0.352.Like the CM, the SLSEM receives higher values of σ′ xy than other methods.The JM correlates best with the SEM, for which the r was 0.711.In turn, the correlation considering the JM, with r = 0.205, was the weakest.On the other hand, the JM correlated best with bi-logarithmic methods.The correlation value between them was determined at the level of 0.719.
To conclude, the results of the correlational analysis revealed that the bi-logarithmic methods correlated best with the other methods of determining σ′ xy .The strongest correlations were observed between the bi-logarithmic methods and the methods based on work approach, i.e.WM and WPUVSM (r = 0.976); slightly worse correlations were archived with the semi-logarithmic methods of BM (r = 0.948) and PSM (r = 0.923) and the DSEM (r = 0.789), SEM (r = 0.762), and JM (r = 0.719).

Agreement between the methods
The difference charts were used to assess the compatibility of particular methods with each other.In the work presented herein, the primary purpose of the analysis was to assess the differences between the results obtained from different methods for determining σ′ xy .According to the division presented earlier, the methods have been assigned to individual groups.Therefore, one method was selected from each group, and its results were the compression behaviour of this kind of soil has been represented by the two straight lines in the plot of ln(1 + e) against log σ′ v .
In the case of interpreting the test results, not all methods were easy to apply, or there were problems with unequivocally determining the reliable value of σ′ xy on their basis.This mainly concerned SEM, SLSEM and JM.www.nature.com/scientificreports/involves considerable judgment during interpretation, a much more straightforward and firmly theoretical approach based on work is recommended.2. Among the other semi-logarithmic methods, the σ′ xy values determined using PSM and BM differed from CM by 65 kPa and 50 kPa, respectively.Thus, the bias between PSM and BM amounted to 15 kPa.3. Considering the bi-logarithmic methods (i.e.BTM, OIM and ONM), complete convergence of the determined σ′ xy values has been observed for all analysed samples.This is due to the same determining procedure and the lack of change in the course of the compressibility curve regardless of the log (1 + e) − log σ′ v or ln(1 + e)-log σ′ v graph used.The same situation pertained to methods based on the concept of work, i.e.WM and WPUVSM.However, it should be borne in mind that methods based on two logarithmic scales may hide the scatter of the recorded data during the test; therefore, it is suggested to concentrate the points on the curve or even smooth its course before applying them.4. The convergent validity analysis demonstrated the highest positive association between bi-logarithmic methods and WPUVSM (i.e.r = 0.976) and bi-logarithmic methods and BM (i.e.r = 0.976), whereby the bias between them was about 36 kPa and 10 kPa, respectively.This means that regardless of the shape of the investigated compression curve, the type of soil or the sampling depth, if one has two related methods, the interpreter who scores, for instance, a high value of σ′ xy by one method should score high by the other as well. 5. Mean difference analysis using difference plots showed that relatively similar σ′ xy values were obtained using the following methods: CM against methods based on work approach (i.e.WM, WPUWSM, and DSEM) and bi-logarithmic methods against BM and JM.In this connection, at least three values of σ′ xy determined by different methods should be averaged and adopted as reliable.Considering the results of the analyses performed, WPUVSM, BM, and BTM can be used.However, it should be remembered that in most comparisons, there were extensive ranges of established agreement limits, which may obscure the unambiguous assessment of all methods.They are all based on more or less complex graphical procedures, and it is practically impossible to say which approach is correct.

Figure 2 .
Figure 2. Location of the sampling area.

Figure 3 .
Figure 3. Classification of the soils used in the present study.

Figure 4 .
Figure 4. Soil compression curves in a semi-logarithmic plane: (a) void ratio vs effective vertical stress, (b) normalised void ratio vs effective vertical stress.

Figure 5 .
Figure 5. Soil compression behaviour: (a) types of compression curve in a normalised void ratio versus effective vertical stress space, (b) compression data plotted in a bi-logarithmic plane.

Figure 8 .
Figure 8. Correlation matrix for convergent validity analysis (Note: *the correlations were performed on a reduced number of samples due to outliers).

Table 1 .
Physical parameters of intact intermediate soils utilised in the present study.F Cl is the clay fraction, F Si is the silt fraction, F Sa is the sand fraction.

Table 2 .
Compression parameters of analysed soils.