Transformation models for the compressibility properties of Finnish clays using a multivariate database

ABSTRACT The compression index obtained from an oedometer test is often used to estimate the settlements of clayey subsoil, but compressibility parameters are rarely available during the preliminary geotechnical design phase. Various empirical correlations linking compressibility to other properties such as water content have been proposed. However, as Scandinavian clays are soft and exhibit greater compressibility, the existing transformation models for compressibility can be biased when applied to Finnish clays. This paper compiles a partial multivariate database of Finnish clayey soils and demonstrates that the existing transformation models tend to underestimate the compressibility of Finnish clays. The new transformation models are constructed by means of a 2-degree polynomial regression applied to the natural logarithms of the soil properties. Finally, the transformation uncertainties are quantified via the standard deviation of errors and the coefficient of variation. The best predictors for the compressibility of Finnish clayey soils were found to be the void ratio and water content. When the void ratio was combined with a secondary predictor, such as the ratio between undrained shear strength and preconsolidation pressure or plastic limit, the transformation uncertainty decreased notably.


Introduction
The compression index C c , which is defined as the decrease in void ratio per tenfold increase in vertical stress, is often used to estimate the compression of clayey subsoil. C c is obtained using an oedometer test performed on an undisturbed specimen. However, as oedometer tests require good-quality specimens and an incremental loading oedometer test usually lasts several days, compressibility parameters are rarely available when the preliminary geotechnical design is being done. Consequently, many researchers have presented empirical regression equations (i.e. transformation models) to estimate C c from other properties such as water content (w n ) or liquid limit (w L ) (Skempton 1944;Helenelund 1951;Hough 1957;Moran and Mueser 1958;Nishida 1956;Cozzolino 1961;Terzaghi and Peck 1967;Elnaggar and Krizek 1970;Azzouz, Krizek, and Corotis 1976;Ogawa and Matsumoto (1978) ;Mayne 1980;Koppula 1981;Saarelainen 1981;Al-Khafaji and Andersland 1992;Kootahi and Moradi 2017;Park and Lee 2011;Di Buò et al. 2019).
Regarding Finnish fine-grained soils, very few transformation models are available for compressibility. Most Finnish geotechnical handbooks present only two models: Helenelund's (1951) non-linear empirical correlation for compressibility index using w n and Janbu's (1988) model, which links the modulus number to w n . However, neither of these transformation models were derived specifically for Finnish clays, and they lack the statistical measures required for use in reliability analyses. Because these transformation models for compressibility are characterised by a rather large transformation uncertainty, it is essential that they can be easily defined in probabilistic terms (e.g. Phoon 2017). Recent studies on the transformation models for compressibility of Finnish clays have been conducted, but they are characterised by a rather small number of distinct sites or a small geographical area. For instance, Saarelainen (1978Saarelainen ( , 1981 created linear regression models for municipal data (Helsinki capital area) and Di Buò et al. (2019) studied the empirical correlations for the compressibility of soft sensitive clays from four sites in Finland.
Updated transformation models for the compressibility of Finnish clayey soils are clearly needed. Thus, the main objectives of this paper are: (i) to compile a large and versatile database of compressibility of Finnish clays and gyttja soils; (ii) to compare the existing transformation models with the compiled database; and (iii) to derive new transformation models for the compressibility of Finnish clayey soils (namely, C c , compression ratio and swelling index).
In this paper a regional database of Finnish clayey soils consisting of various geotechnical properties is presented. The compiled database contains 856 oedometer test results from over 30 sites. The database is partially multivariate, including up to 14 parameters of interest; for instance, almost half of the oedometer tests were accompanied by the undrained shear strength of an adjacent specimen. The database is called FI-CLAY/ 14/856, following the name format recommended by TC304 (http://140.112.12.21/issmge/tc304.htm) and Phoon (2019). To the authors' knowledge, this is the largest database of oedometer tests on Finnish clays to date. Another extensive multivariate database of Finnish clays was compiled by D' Ignazio et al. (2016), but it does not include compressibility properties. One practical benefit of multivariate databases is that they can be used to derive multivariate distributions. Effectively, the uncertainty in one geotechnical parameter is reduced when one or more additional parameters are available (Müller, Larsson, and Spross 2014;Ching and Phoon 2014a;Phoon et al. 2016).
This paper is organised as follows. Firstly, a brief overview of the existing transformation models for compressibility is presented. Secondly, the complied database of Finnish clays FI-CLAY/14/856 is described and compared to existing transformation models. Next, any potential outliers are removed and new transformation models for compressibility are derived for FI-CLAY/14/822 using polynomial regression applied to the natural logarithms of the soil properties. Finally, the transformation uncertainty is evaluated using statistical measures.
2. Overview of the existing transformation models for compressibility 2.1. Measures of compressibility C c is commonly used to calculate one-dimensional vertical strain, which is defined as the change in soil layer thickness divided by the layer thickness. In the case of a normally consolidated (NC) clay layer, in which preconsolidation pressure σ p ' equals effective in-situ stress σ v0 ' , the strain ε NC is calculated using Equation (1): where Δσ v is the added vertical load and e 0 is the initial void ratio. The term C c /(1+e 0 ), known as the compression ratio (CR), is considered by some researchers to be a better measure of compressibility than C c (e.g. Janbu 1988;Wesley 1988). Similarly, for over-consolidated (OC) clays (i.e. σ p ' greater than σ v0 ' ), the strain ε OC+NC is calculated using Equation (2): where C s is the swelling index, which is defined as the slope of the swelling line (or recompression line). By dividing C c and C s by ln(10) (≈2.3), one acquires approximations of the Cam-Clay parameters λ and κ, respectively (ln[10] accounts for conversion to natural logarithm). The typical shape of oedometer curve for soft Finnish clays together with the interpretations for both C s and C c is shown in Figure 1. The stress-strain behaviour of many marine Finnish clays is so non-linear in semi-logarithmic scales that the C c method cannot be applied without precautions (e.g. Tanaka 2002;Di Buò et al. 2019). The common practice is to determine C c using two load increments after the preconsolidation pressure, which usually corresponds to the relevant stresses in a practical geotechnical design. Due to the potential limitations of C c , the Ohde-Janbu method is often used in Finland instead. In this model, the strain ε NC caused by the stress increase from preconsolidation stress σ' p to final stress (σ v0 ' + Δσ v ) was calculated using Equation (3) (Janbu 1963;Helenelund 1974): In Equation (3), modulus number m 1 and stress exponent β 1 are the parameters defined from oedometer tests by means of curve fitting and σ ref is the reference Figure 1. The interpretations of C c and C s for soft clays. stress (100 kPa). If β 1 = 0, the stress-strain curve is linear in a semi-logarithmic scale, as in the C c concept. As a result, C c equals (1+e 0 )ln(10) divided by m 1 . Equation (3) can also be used to calculate the strain in an overconsolidated state ε OC (stress states from σ vo ' to σ' p ), in which case the parameters are usually marked with subscript "2" (m 2 and β 2 are defined using the unloading-or reloading-curve of the oedometer test). Although Ohde-Janbu is usually more accurate for modelling the compressibility of soft Finnish clays (due to closer replication of the oedometer curve, see Figure 1), its parameters cannot be used as the target parameters of simple transformation models due to their correlation; the Pearson's coefficient of correlation is ρ = 0.40 for m 1 and β 1 (see Table S2 in the online supplementary material).

Existing transformation models for compressibility
Most of the existing transformation models for compressibility of clays are simple linear regression models with a single predictor. The earliest models for C c were constructed using w L as the predictor (Skempton 1944). However, subsequent studies showed that using the initial void ratio e 0 as the predictor leads to a higher R 2 value (coefficient of determination) (e.g. Azzouz, Krizek, and Corotis 1976;Al-Khafaji and Andersland 1992;Kootahi and Moradi 2017). Similarly, w n is associated with high R 2 , although Al-Khafaji and Andersland (1992) pointed out that many studies have reported w n values defined for partially saturated soil specimens (rather than fully saturated), which in turn can lead to weaker model accuracy. In Finland, the specimens used in an oedometer test are almost always fully saturated; consequently, w n is often selected as the predictor for C c (Helenelund 1951;Saarelainen 1981Saarelainen , 1978Di Buò et al. 2019). Indeed, Di Buò et al. (2019) observed that the relationship between C c and w L , plasticity index PI, or clay content is characterised by larger scatter compared to water content w n . Likewise, Kootahi and Moradi 2017 argued that single-predictor models incorporating index properties (i.e. w L or PI) are characterised by larger scatter due to greater measurement uncertainty compared to physical properties such as w n or e 0 .
Many researchers recommend using a two-predictor model with void ratio e 0 (or w n ) combined with w L (or PI) to account for both the structure and mineralogical effects on clay compressibility (e.g. Al-Khafaji and Andersland 1992; Kootahi and Moradi 2017). That is, the mineralogical composition and particle-surface characteristics can be represented by intrinsic index properties such as w L or PI, whereas the physical properties (e 0 and w n ) represent the structure induced by the sedimentation environment. Hence, depending on the clay deposit characteristics, adding w L or PI as an additional predictor may increase the model accuracy. On the other hand, Azzouz, Krizek, and Corotis (1976) concluded that including w L as an additional predictor did not significantly increase the accuracy of C c or CR predictions.
Besides w L or PI, other additional predictors has been suggested as well. For instance, Leroueil, Tavenas, and Bihan (1983) suggested adding sensitivity S t to the relation between e 0 and C c . Al- Khafaji and Andersland (1992) argued that the multiple parameter models are more accurate especially for organic clays, as earlier studies (e.g. Al-Khafaji 1979) have showed that the C c of organic soils is highly dependent on organic content and preconsolidation pressure. Some existing transformation models for compressibility (i.e. C c , CR and C s ) are presented later in Table 2. 3. Compiled clay database (FI-CLAY/14/856)

Collected multivariate data
The database of Finnish clays (FI-CLAY/14/856) compiled for this study contains mostly glaciomarine finegrained soils with some specimens from glaciolacustrine origins. The full database can be found in the online supplementary material (Table S1). The compiled database was based on two sources: Subset A (Table A1 in the Appendix), which consists of multivariate data from literature, and Subset B (Table A2), which includes data collected for this study specifically (35% of all data). Subset A consists of three sites from Gardemeister (1973) (3% of all data), data compiled by the City of Vantaa (Pajunen 1976) (13% of all data), and the RITA database (Rekonen and Lojander 1999) (49% of all data). The RITA database consists of oedometer tests performed and analysed by the geotechnical research group at the Helsinki University of Technology (now Aalto University).
Most sites in the compiled database FI-CLAY/14/856 represent small areas or test embankments; in other words, the distance between sampling profiles is in the scale of 3-30 m. At the Suurpelto (Subset B) and Vantaa sites (Subset A; Pajunen 1976) however, the maximum distance between the sampling profiles was approximately 1-3 kilometres. The Söderkulla-Nikkilä (RITA database) specimens were collected as a part of a road renovation project, so their maximum separation distance is up to 10 kilometres.
Basic statistics for the parameters of interest (14 soil properties) are presented in Table 1. Each observation (table row) in the FI-CLAY/14/856 database represents soil properties at a specific depth and location. Each data point (n = 856) has at least three soil properties: w n , e 0 and C c . Other data points have additional soil properties to a varying extent; in other words, the database is partly multivariate. For example, Table 1 shows that 413 data points contained both the undrained shear strength s u (determined from a fall cone test) and C c (i.e. the oedometer test result). Table 1 also shows the corresponding soil property statistics for a global database of marine fine-grained soils compiled by Kootahi and Moradi (2017). These databases are also compared in Figure 6. Table 1 shows that the main differences between the regional and global databases are the following: the global database contains higher values of w n , e 0 and PI whereas the FI-CLAY/14/856 contains higher values of C c , CR and S t . These Finnish sensitive clays may be one explanation for higher values of C c (e.g. Leroueil, Tavenas, and Bihan 1983;Di Buò et al. 2019).
Data from Gardemeister (1973) and Pajunen (1976) were originally fully multivariate, and the RITA database was partially multivariate. For the remaining data, multivariate data points were interpreted for this study. The oedometer test results and other properties (classification and strength properties) were combined if both of the following applied: (i) the vertical distance between specimens taken from the same sampling point was not more than 100 mm; and (ii) the w n of the oedometer test specimen and the other specimen did not differ more than 10 percentage points. It should be noted that as some of the oedometer tests were performed on specimens within a 100 mm section, there are a few multivariate rows which include duplicate values of the other property (e.g. plastic limit).

Soil types
In the FI-CLAY/14/856 database, the specimens were classified according to the Finnish geotechnical classification system (Korhonen, Gardemeister, and Tammirinne 1974). In this system, the grain size distribution and organic content defined the soil name, meaning that plasticity was not considered and thus many of the data points lack consistency limits. The naming convention based on clay fraction (grain size smaller than 0.002 mm) and organic content are presented in Figure 2 (a), together with the number of each soil type in the database. The pie chart in Figure 2(a) shows that most of the specimens were clays. For some specimens (i.e. "fine-grained soils"), the soil type was not reported.
The specimens in the database which were accompanied with consistency limits and complete soil classification are plotted in the plasticity chart in Figure 2(b). According to the Unified Soil Classification System (ASTM D 2487-06), most of the specimens (Fat clays, Organic clays, Organic silts, and Gyttja) lie above the A-line, i.e. they represent CH or OH types. On the  Kootahi and Moradi (2017) for model building. Notes: n = number of data points; COV = coefficient of variation; w n = natural water content (%); e 0 = initial void ratio; w L = liquid limit (%); F = fineness number (%); w P = plastic limit (%); γ = unit weight (kN/m 3 ); Org = organic content (%); Cl = clay fraction (%); s u = undrained shear strength (kPa) (fall cone); S t = sensitivity; σ' p = preconsolidation pressure (kPa); OCR = over-consolidation ratio; C c = compression index; CR = compression ratio; C s = swelling index.
other hand, the majority of Lean clays and Clayey silts represent CL or OL types, and a few Clayey silts fall into the CL-ML type. Finally, Figure 2(c) represents the C c (w n ) relationship for different soil types. Gyttja specimens (including Clayey Gyttja) behaved somewhat differently compared to other soil types, but there was also overlap with other soil types (i.e. it was not necessary to remove the Gyttja specimens from the data).

Oedometer test types
Most of the data (n = 723) represent incrementally loaded (IL) standard oedometer tests with a 24-hour load sequence and loading which is doubled for each increment. However, the data also include some speciality oedometer test types, which are shown in the C c (w n ) graph presented in Figure 3: "horizontal" tests were performed on specimens rotated 90 degrees, "permeability" IL tests included a falling head test performed at selected load increments and "creep" tests contained longer-than-usual load increments. However, for most of these creep tests, the strains at 24 h were utilised instead of the final values. Besides 24-hour values, most of the recent oedometer test curves were constructed using the end-of-primary (EOP) values; those of the normal compression line ( Figure 1) were estimated by means of consolidation time analysis (e.g. the square root time method) in order to minimise the effects of creep on the interpreted parameters. Figure 4 illustrates these EOP results (n = 128) in a C c (w n ) graph. All oedometer tests in Subset B sites except the Otaniemi test site and the POKO site represented EOP results. Figure 3 also illustrates the constant rate of strain (CRS) oedometer tests (n = 133). Three sites contained C c(+10kPa) = 0.000102(w n ) 2.284 856 0.54 0.51 C c(+20kPa) = 0.000289(w n ) 2.031 856 0.57 0.46 C c(+50kPa) = 0.0011(w n ) 1.669 856 0.73 0.43 Moran and Mueser (1958) Organic soils, peats, organic silts and clays C c = 0.0115w n 856 1.34 0.49 Kootahi and Moradi (2017) Marine fine-grained soils  Kulhawy and Mayne (1990) All clays C s = 0.0027PI 267 1.18 0.87 Tiwari and Ajmera (2011) Reconstituted montmorillonite dominated soils C s = 0.0087(e 0 ) 2 + 0.031(e 0 ) 656 0.96 0.43 C s = 6 · 10 −6 (PI) 2 + 0.0026(PI) 267 1.15 0.93 Notes: n = number or samples; b = mean bias; δ = COV of the bias values; C c = compression index; CR = compression ratio C c /(1 + e 0 ); C s = Swelling index; w n = natural water content (%); e 0 = initial void ratio; w L = liquid limit (%); PI = plasticity index.
CRS tests: Haarajoki (n = 34), Söderkulla-Nikkilä (n = 52) and Perniö test embankment (n = 47). The applied strain rate varied greatly due to rate-dependency studies (0.03-15%/hour). The OCR values of CRS test specimens were calculated using the apparent preconsolidation pressure (i.e. no rate correction was applied). Figures 3 and 4 demonstrate that none of these non-standard oedometer tests formed their own distinct populations and therefore they may be used to derive the transformation models for compressibility. It should be noted that Figure  4 also illustrates specimens with low OCR which are further discussed in "Removal of outliers" section.

Applied approximations
As explained above, in Finland the Ohde-Janbu method is generally preferred over C c , so many of the oedometer tests lack interpretations of C c : For these tests C c was estimated using the Ohde-Janbu parameters m 1 and  β 1 . This approximation was performed using the method proposed by Löfman and Korkiala-Tanttu (2019) and applied to the Subset A data and all the CRS data from the RITA database ("C c estimated" in Figure 4). In Gardemeister's (1973) data, e 0 was approximated by assuming full saturation and density of solid particles ρ s = 2.7 g/cm 3 (i.e. e 0 = 2.7(w n /100)). Because these specimens were taken from depths below the groundwater level, the estimation based on w n is quite accurate. Figure  S1 in the online supplementary material shows that the relationship between w n and e 0 is linear with little scatter. Accordingly, the Pearson's coefficient of correlation for w n and e 0 was found to be very high (ρ = 0.98) (see Table S2 in the online supplementary material).
Finally, some data points contained fineness number F. This parameter was determined from the fall cone test on remoldedspecimens, and the result was taken as the approximation of w L determined using the Casagrande test (e.g. Gardemeister 1975). For many specimens, both parameters were defined, and it was verified that their relationship was linear and almost non-biased (see Figure S2 in the online supplementary material). Likewise, the ρ for w L and F was found to be as high as 0.98 (see Table S2 in the online supplementary material).

Comparison of existing transformation models to FI-CLAY/14/856
In this section, the database is compared to existing transformation models for compressibility. The bias factors and transformation uncertainties are evaluated for the selected existing models. The bias b i of ith data point can be defined as the actual target value (C c , CR, or C s ) divided by the predicted target value. The bias factor b is the sample mean of all biases b i . Thus, if b = 1, the model is unbiased. The sample coefficient of variation of values b i is the COV value of the transformation model, denoted by δ. According to of Ching and Phoon (2014b), the variability term ε for a transformation model can calculated using Equation (4): The mean of ε is 1 and its sample coefficient of variation describes the deviation between the actual target value and the unbiased prediction. The COV of ε corresponds to δ. If δ = 0, there is zero data scatter about the transformation model, i.e. there is zero transformation uncertainty.
The bias factors b and δ (COV of ε) for the selected existing models are collected to Table 2. The selected models are also illustrated in Figures 5-9 together with the database values. Transformation models which have been found to be inapplicable to marine fine-grained soils (Kootahi and Moradi 2017) have been omitted from the calibration. In addition, predicted target values, which were either negative or smaller than the observed minimum value, were removed from calibration data. Calibration results in Table 2 show that most of the existing transformation models compressibility had bias factors greater than 1, meaning that the target value was generally underestimated.  However, models with multiple predictors had smaller biases on average. Next, each major predictor for compressibility is discussed in detail. Figure 5 illustrates that the existing C c (w n )-transformation models performed well for soil specimens with low w n (less than 75%), whereas all the existing modelsexcept Di Buò et al.'s (2019) models for Finnish sensitive claysunderestimate C c for specimens with higher w n . Furthermore, the upper range reported by Helenelund (1951) overlapped the data. The overestimation by Di Buò et al.'s (2019) models is due to different interpretation method applied (C c interpreted from small stress increments after the preconsolidation pressure in a CRS test) and due to high sensitivity of the studied clay specimens. In fact, as discussed later, some of the high C c values in FI-CLAY/14/856 represent sensitive Perniö clay included in both of these datasets (i.e. different tests results for the same clay layer).
A similar under-prediction of C c was observed using transformation models with e 0 ( Figure 6) and w L (Figure 7). Figure 7 also includes the above-mentioned fineness number F, which can be considered as the approximation of w L . Moreover, the relationship between C c and w L is characterised by significant scatter, which also explains the high transformation Figure 7. Existing C c (w L ) models. Figure 9. Existing C s (w n ) models and the newly derived model.  uncertainties (δ in Table 2). Table 2 shows that neither w L nor PI should be used as the single predictor of compressibility for FI-CLAY/14/856 due to high values of δ. Finally, Figure 8 shows the data together with the selected transformation models for CR. As observed above, the existing models performed well at lower w n values but their biases increased with w n .
The existing models for estimating C s with w n are illustrated in Figure 9, which shows that the existing model for Finnish sensitive clays provided by Di Buò et al. (2019) overestimated C s . The reason for this overestimation could be the testing conditions, as Di Buò et al. (2019) actually defined the initial recompression index (prior surpassing the preconsolidation pressure) rather than C s at unloading-reloading. Moreover, they observed that the coefficient of determination between w n and the recompression index was very weak (R 2 = 0.16). Similarly, the C s data (FI-CLAY/14/856) is marked with greater scatter compared to C c and CR. Nevertheless, the model suggested by Kootahi (2017) fits rather well with the data, and the bias factor b is approximately 1 ( Table 2). On contrary, the relationship between PI and C s was characterised by larger scatter (see Figure S3 in the supplement); accordingly, calibration yielded very high δ values.

The polynomial regression model
The calibration studies above demonstrated that most of the existing transformation models underestimate the compressibility of Finnish clays. Moreover, as the scatter increases with predictor value (e.g. Figure 5), a simple linear regression model cannot provide a robust definition of the transformation uncertainty. Hence, a 2degree polynomial regression model applied to the natural logarithms of the soil properties (i.e. the target and the predictors) was selected. The polynomial model, which is an extension of the linear regression, was fitted using the ordinary least squares (OLS) method within the Scikit-learn library for Python (Pedregosa et al. 2011). The transformation model (a paraboloid) was defined using Equation (5): whereŶ is the predicted ln-transformed target (such as ln[C c ]), X 1 is the ln-transformed primary predictor (such as ln[e 0 ]) and X 2 is the ln-transformed secondary predictor. The b terms are regression coefficients. In the case of single-predictor model, the coefficients for variables containing X 2 would be zero. The reason for using a natural logarithm (ln-transformation) is explained below. The transformation error ε is a zero-mean normal random variable with standard deviation σ ε . For a single data point, this transformation error ε i was quantified as the residual (error) term of ith observation, i.e. the actual value Y i defined from the oedometer test result minus the predicted valueŶ i . The standard deviation of the transformation error σ ε can then be estimated from these individual residual errors using Equation (6) (e.g. Ang and Tang 2006): where n is the number of samples used in the linear regression and v is the number of degrees of freedom (that is, the number of estimated regression coefficients).
The transformation model presented above was selected in order to satisfy two assumptions usually made in OLS linear regressions: (1) Homoscedasticity, which means that the standard deviation of the transformation error σ ε (Equation (6)) should be equal across the regression line. After applying the ln-transformations and fitting the polynomial regression model, the scatter became almost constant along the predictor's value (e.g. Figure 12(a)).
(2) The normality of the residuals ε i . After applying the abovementioned transformations and polynomial fit, the residuals became normally distributed (see Figure S4 in the online supplement).
Another assumption in OLS to be considered is the independence of the predictors: If the predictors are strongly correlated, a phenomenon known as multicollinearity can lead to unstable regression coefficient estimates. Predictors such as w n and e 0 were thus not combined as predictors in this paper.
The standard deviation of error σ ε provides one measure of the transformation model uncertainty but it cannot be used to compare models for different target parameters. Hence, the definition by Ching and Phoon (2014b) was used to derive another measure of transformation uncertainty, i.e. the coefficient of variation δ of the transformation error (Equation 4). Hence, the value of δ was derived similarly as in the calibration analysis. Even though δ was originally intended for calibration studies, its value is an indication of the transformation model quality (Phoon 2017). In addition, the mean and median of bias values b i were calculated for each model to confirm that the newly derived models are unbiased.

Removal of outliers
Potential outliers have to be detected and removed for two reasons: (1) they may distort the fitting of the transformation models, thus resulting in a biased model; and (2) the estimated transformation uncertainty should include a minimal number of other sources of uncertainty, such as measurement error or inherent variability (Phoon and Kulhawy 1999). Outliers related to measurement error can result from factors such as corrupted data (e.g. human error in data transfer) or specimen disturbance. On the other hand, the uncertainty related to inherent variability is caused by the very nature of multivariate databases; only a limited number of laboratory tests can be performed on a single soil specimen, and no two soil specimens are perfectly identical. Finally, the data point may come from another population entirely; for example, soils with very high organic content might behave like outliers due to their extreme physical properties.
First, under-consolidated and potentially disturbed specimens are detected using the physical constraints of OCR (over-consolidation ratio, σ' p /σ' v0 ). Specifically, data points with an OCR smaller than 0.7 were removed as outliers (29 tests were conducted in total; see Figure 4). The threshold value of 0.7 was selected because the OCR for normally consolidated Finnish clays varies between 0.8 and 1.5 (Korhonen, Gardemeister, and Tammirinne 1974). Finally, the threshold value of 0.7 (instead of 0.8) was selected because the value of OCR is also affected by uncertainty in determining the preconsolidation pressure from incrementally loaded odometer tests.
Next, the extreme values of C c were identified. A histogram of C c (Figure 10) showed that there were five oedometer tests which produced very high C c . In particular, with a bin width of 0.5, all other bins (C c < 5) contained at least two records. Based on this observation, tests with C c ≥ 5 (five tests in total) were considered rare cases and thus outliers. In fact, these tests were all CRS tests on highly sensitive Perniö clay (with the same sampling point and the same depth) and the C c slope was interpreted as the tangent positioned right after the preconsolidation pressure, resulting in an almost vertical slope. To summarise, in total 34 data points (4% of all data) were considered outliers and removed prior to deriving the transformation models for C c and CR. In the case of C s , a similar extreme-value criterion was applied and an additional three records (C s > 0.4) were removed as outliers ( Figure 10).
Finally, another outlier detection was applied during each transformation model fitting separately to account for factors such measurement error and inherent variability. This detection was based on the standard deviation method with a 3σ-threshold. In other words, those data points in which the residual error ε i was more than three standard deviations (σ ε ) apart from the zero-mean, were detected as outliers and hence removed before defining the final regression coefficients. The number of outliers removed in this last stage are presented in Tables 3 and 4 and in Figure 12(a), Figures 13-16(b). In addition, the figures illustrate the fitted model and adjusted R-squared (R 2 adj ) also for the datasets including these outliers.

Leave-one-group-out analysis
Because some of the sites in the database contained significantly greater numbers of data points than others, a leave-one-group-out (LOGO) analysis was performed in order to verify that no single site had too great an influence on the transformation model. A simple single-variable linear model for ln(C c ) with ln(w n ) was selected for the LOGO study. In the LOGO analysis, one group (here one site) is left out of the training data each time; that is, the linear regression coefficients are estimated (fitted) using all the remaining data except that one group. The LOGO analysis recognised three sites (Söderkulla-Nikkilä, Haarajoki and Suurpelto) Figure 10. Histograms of C c and C s . that caused the most significant changes to the regression coefficients. C c (w n )-transformation models (polynomial regression, Equation (5)) were then derived for all data and for the datasets without each of these extreme cases. Figure 11 shows that the difference in the models is hardly noticeable. At higher values of w n , the difference in the models was greater, but so was the scatter in C c . This indicated that all the sites can be included in the fitting of the transformation models without the risk of creating a biased model due to a large sample from another population.

New transformation models and their transformation uncertainty
The most suitable predictors for the transformation models for C c and CR were identified using a correlation matrix (see Table S2 in the online supplementary material). The three most suitable predictors were found to be w n , e 0 and unit weight γ. The estimated regression coefficients and the model metric for these single-predictor models are presented in Table  3. The corresponding scatter plots for models with w n are presented in Figure 12(a,b) and Figure 13(a,  b). The fitted models for CR and C c are also illustrated in Figures 8 and 11, respectively. All the other scatter plots are presented in Figures S6-S7 in the online supplementary material. For both C c and CR, the models with e 0 had the smallest transformation uncertainty; for C c (e 0 ) model, the R 2 adj was 0.83, the standard deviation of error σ ε was 0.346, and the coefficient of variation δ was 0.35. Using w n as the predictor resulted in very similar model metrics. On the other hand, the models with w L and PI as single predictors resulted in clearly higher transformation errors, and hence these results have been omitted from Table 3. For example, for C c (w L ) model, the R 2 adj was only 0.38. Table 4 presents the regression coefficients and model metrics for two-predictor models for C c and CR. Void ratio e 0 was selected to be the primary predictor due to the lowest transformation uncertainty observed amongst single-predictor models. The four most suitable secondary predictors were found by studying the correlation between the potential secondary predictors X 2 and the errors ε i of the C c (e 0 ) model. A high correlation would effectively imply that this another predictor could increase the accuracy of the model if added as a secondary predictor (Ching and Phoon 2014b). The highest ρ values were found with ratio (s u /σ p ' ) (ρ = −0.29), preconsolidation pressure σ p ' (ρ = 0.20), plastic limit w P (ρ = −0.21), and organic content Org (ρ = −0.19). On the other hand, the corresponding correlations with w L and PI were ρ = −0.03 and ρ = 0.01, respectively. Similar correlation results were acquired for the errors of the CR(e 0 ) model. Models with three predictors were also derived, but the decrease in transformation uncertainty was considered to be too Table 3. Single-predictor transformation models for compressibility of Finnish clays.
Model variables Outliers Regression model: 38 Notes: C c = compression index; CR = compression ratio C c /(1 + e 0 ); C S = Swelling index; n = number of samples; R 2 adj = adjusted R-squared; σ ε = standard deviation of the transformation error; MD = median; δ = COV of the bias values; w n = natural water content (%); e 0 = initial void ratio; γ = unit weight (kN/m 3 ). Figure 11. New C c (w n ) model and the models fitted to data without three most influential sites.
small to outweigh the potential uncertainties brought by added complexity of the model.
Graphs comparing the predicted target values and the actual values are illustrated for one secondary predictor: C c (e 0 , w P )-model is presented in Figure 14 and CR(e 0 , w P ) in Figure 15. Adding these secondary parameters resulted in a notable decrease in transformation uncertainty. For example, for C c (e 0 , w P ) model, the adjusted R 2 adj was 0.88, σ ε was 0.317, and δ was 0.32. (For the same data, the metrics for the C c (e 0 ) model were R 2 adj = 0.85, σ ε = 0.350, and δ = 0.36.) The smallest transformation uncertainty was reached by using the ratio (s u /σ p ' ) as the secondary predictor. On the other hand, using w L as the secondary predictor improved the model metrics only slightly (R 2 adj = 0.85, σ ε = 0.350, and δ = 0.36), and hence this model was omitted from Table 4.
For swelling index C s , the transformation uncertainties regarding the single-predictor models were greater compared to models for C c and CR (Table 3). Overall, R 2 adj was less than 0.50 for C s models. The scatter plots for C s (w n ) model are shown in Figure 16(a,b), and the  fitted model is also illustrated in Figure 9. The remaining scatter plots can be found from the online supplement ( Figure S8). Regarding the potential secondary predictors for the C s (e 0 ) model, the highest ρ value was found for clay content (ρ = 0.14). However, adding this secondary predictor did not notably increase the model metric, and hence two-predictor models for C s are omitted from Table 4.
Considering the results discussed above, the following equations can be recommended for estimating compressibility of Finnish fine-grained soils:

Discussion
Contrary to previous findings on multiple predictor models for clay compressibility (e.g. Kootahi and Moradi 2017), adding w L or PI as the secondary predictor did not notably decrease the transformation uncertainty in estimating C c of Finnish fine-grained soils. As suspected by Di Buò et al. (2019), the physical properties (such as w n ) seem to better explain the variability of C c for Finnish clays. On the other hand, adding the plastic limit w P as the secondary predictor did result in a notable decrease in transformation uncertainty for C c .  34 Notes: C c = compression index; CR = compression ratio C c /(1+e 0 ); n = number of samples; R 2 adj = adjusted R-squared; σ ε = standard deviation of the transformation error; MD = median; δ = COV of the bias values; e 0 = initial void ratio; σ' p = preconsolidation pressure (kPa); w P = plastic limit (%); s u = undrained shear strength (kPa) (fall cone); Org = organic content (%).  As such, for Finnish clays, w P may be a better representation of the non-structural effects (such as mineralogy and particle-surface characteristics) to compressibility. Kootahi and Moradi (2017) also used w P as the secondary predictor of C c for marine fine-grained soils, but the increase in model metrics was greater when w L or PI was used instead.
The smallest transformation uncertainty was reached by using the ratio (s u /σ p ' ) or σ p ' as the secondary predictor. In a preliminary design stage, σ p ' is usually not available, and hence these transformation models are intended to be used in Bayesian approaches. For instance, one site-specific oedometer test is enough to compute the prior distribution of C c for Bayesian updating via the reported standard deviation of error σ ε . Alternatively, the engineer can utilise the limited sitespecific data to verify whether the predicted C c is within the 95% confidence limits as instructed by Ching (2017). In general, the transformation uncertainty was found to be moderate at δ = 0.30-0.43. On the other hand, transformation models for the undrained shear strength of clays exhibit somewhat smaller transformation uncertainty (e.g. δ = 0.25-0.30 for the most reliable models derived by D'Ignazio et al. [2016]). Higher transformation uncertainty for compressibility is to be expected, since the interpretation of C c from the oedometer curve of soft clays is not explicit. In practice, the 95% confidence interval for a settlement prediction can be approximated by considering the ±2σ ε limits of the predicted C c ; for a five metres thick clay layer with void ratio of e 0 = 1, a two-metre high embankment (corresponding to 40 kPa) on normally consolidated clay would result in a settlement of 77-308 mm (assuming effective in-situ stress of 40 kPa). The prediction's reliability decreases with increasing values of e 0 and w n . On the other hand, a ground-supported foundation is not usually considered a feasible option for highly compressive soft clays, hence reducing the need for accurate settlement predictions.
The transformation uncertainty related to swelling index C s was notably greater compared to C c and CR. The greater scatter about the trend could be explained by the stress-dependency of C s ; since the unloadingreloading loop of Finnish clays is often marked with hysteresis, the value of C s depends on the stress path applied in the oedometer test (e.g. Vipulanandan et al. 2008).
The statistics of the new models' bias values (Tables 3 and 4) show that the mean bias is systematically slightly larger than 1 (1.04-1.09). Median, on the other hand, is close to 1. This pattern could be explained by the fact that the bias (i.e. the actual target value divided by the predicted value) is log-normally distributed. In fact, the transformation error is often assumed to be lognormal (e.g. Ching 2017). Indeed, the biases for C c (w n ) model seemed to be log-normally distributed based on the histogram's shape (see Figure S4 in the online supplementary material). Other possible explanations for the average bias other than 1 are inaccuracies related to the transformation model or factors such as outliers which remained undetected.
Lastly, it can be concluded that the selected 3σ threshold is a rather strict outlier criteria for geotechnical data, because inherent variability and measurement error also contribute to the observed transformation uncertainty. However, since the standard deviation method is not the most robust outlier detection method (because the statistics which define the outlier criteria are calculated with the very same data which contains the potential outliers), a stricter threshold was selected for this study. Nevertheless, it can be concluded that the apparent transformation uncertainty would be smaller if one would use a less strict (e.g. 2σ threshold) outlier criteria instead.
Finally, it should be noted that the recommended new transformation models (Equations [7]-[11]) should not be used for soil specimens with very high water content (w n > 150%). Figure 8 illustrates how the trend curve for the CR(w n ) model plateaus at highest w n values. Similarly, Figure 13(b) depicts how the predicted CR has a maximum value. Hence, at very high w n values, the transformation model becomes biased. However, as previously discussed, the accuracy of the model is more essential for clays with moderate to low compressibility (i.e. clays with lower w n ).

Conclusions
The following conclusions can be drawn: . An extensive, partly multivariate clay database FI-CLAY/14/856 was compiled and compared to existing transformation models for compressibility of clays. . New transformation models for compressibility of Finnish clays were derived via 2-degree polynomial regression and their transformation uncertainty was evaluated. . The best predictors of compressibility were void ratio and water content. For compression index, the adjusted R-squared was 0.82-0.83, whereas models for compression ratio and swelling index had lower R-squared values (0.67-0.69 and 0.48-0.49, respectively). . When the void ratio was combined with a secondary predictor, such as the ratio between undrained shear strength and preconsolidation pressure or plastic limit, the transformation uncertainty decreased notably for compression index and compression ratio. For compression index, the adjusted R-squared increased to 0.88-0.89. For swelling index, adding a secondary predictor did not notably decrease the transformation uncertainty.