Original Article

Assessing the viscous compression parameter for land subsidence modelling from material characteristics

Pepijn van Elderen^1,2, Gilles Erkens^1,2,3, Cor Zwanenburg^2,4, Harry van Essen² and Esther Stouthamer¹

¹Department of Physical Geography, Utrecht University, Utrecht, The Netherlands; ²Deltares Research Institute, Utrecht, The Netherlands; ³Ecology and Biodiversity, Institute of Environmental Biology, Utrecht University, Utrecht, The Netherlands; ⁴Department of Geoscience & Engineering, TU Delft, Delft, The Netherlands

Abstract

Compaction of soft subsurface layers occurs through different mechanisms, namely consolidation, viscous compression, autocompaction, and shrinkage. Viscous compression is a key long-term process often the most significant component of compaction over long time periods. The rate of viscous compression over time is defined by the viscous compression parameter (C_α), which describes the amount of strain per logarithmic unit of time with constant effective stress. Overall compressibility of different material can often be linked to their material characteristics, yet for viscous compression this is still problematic. Understanding the relationship between material characteristics and C_α is essential for improving compaction predictions in soft subsurface layers such as clay and peat. This study investigates how C_α varies with subsurface material characteristics and explores its predictability based on specific material characteristic values.

An extensive dataset of one-dimensional compression tests on Holocene peat and clay samples from the Rhine-Meuse delta and coastal zone in the Netherlands was analysed. Using principal component and clustering analyses, four distinct lithological classes were identified and labelled based on organic fraction: Clay, Clayey peat, fibrous Peat & decomposed Peat. All four lithological classes show increasing median sample C_α values with decreasing bulk densities and increasing water content and void ratio, with C_α values ranging from 0.0074 in the ‘Clay’ material type to 0.0293–0.0384 in the two ‘Peat’ material classes. Compared to currently applied traditional lithological classifications and default model values in the Netherlands, the material type-specific C_α values provide a more accurate and data-driven alternative. Regression analysis using the full dataset, with water content and dry bulk density as the predictors, yielded a better prediction of C_α (R² = 0.61) than for regressions per material type. Due to large variability within the peaty material types, upper and lower C_α bounds (0.02–0.06) are recommended in practical applications for ‘Clayey Peat’ and ‘Peat’ subsurface layers, as these values are consistently observed across the full range of water content and dry bulk density within these clusters. While the C_α values for ‘Clay’ are similar to standard values (0.0035–0.0123), C_α values for ‘Clayey Peat’ and ‘Peat’ are up to three times higher than the standard values, which would result in twice or thrice the amount of viscous compression over the same time period. Finally, we emphasise the need for adaptable C_α values in models to reflect evolving material characteristics over time, driven by decomposition.

Keywords: classification; engineering geology; Holocene; parameterisation; surficial geology

 

 

Introduction

Land subsidence is often caused by the compaction of subsurface material after fluid extraction, excessive water drainage, loading from constructions, or natural processes and occurs through physical, chemical and biological processes (e.g. Minderhoud et al., 2017; van Asselen et al., 2009; van Elderen et al., 2025). The compaction of subsurface material involves three primary mechanisms: stress-dependent elastic deformation, stress-dependent plastic deformation and time-dependent deformation (e.g. Mesri, 1973; Mitchell & Soga, 2005; Terzaghi, 1941; van Asselen et al., 2009). Time-dependent deformation, referred to as viscous compression, creep or secondary compression, unfolds under sustained loading conditions (e.g. Bjerrum, 1967; Chai et al., 2012; Den Haan & Edil, 1994; Le et al., 2012; Zhao et al., 2020). This slow, continuous compaction, called viscous compression from here onwards, has a particularly large contribution to compaction in clay-rich and organic subsurface layers, which have a high initial compressibility and porosity (e.g. Den Haan & Edil, 1994; Edil & Dhowian, 1979; Le et al., 2012; Mesri & Ajlouni, 2007). Viscous deformation originates from structural adjustments within the matrix of the material that occur over timescales from hours up to decades because of water diffusion from micropores, changes in the adsorbed water layer, and shifting particle interactions (Le et al., 2012; van Elderen et al., 2025). The relative contribution of these driving mechanisms of viscous compression depends on material characteristics and environmental conditions, as well as on decomposition of organic matter, the latter being an important control on the viscous compression mechanisms and structural integrity of peat (Augustesen et al., 2004; Le et al., 2012; O’Kelly & Pichan, 2013; van Elderen et al., 2025).

Numerical modelling plays a pivotal role in quantifying and predicting land subsidence resulting from compaction of unconsolidated subsurface material, i.e. material which has not yet lithified. Existing conceptual models are based on different representations of compaction mechanisms in such materials, such as the nucleus of strain model (Geertsma, 1973; Fokker & Orlic, 2006), the NATSUB3D model (Xotta et al., 2022), constitutive numerical models based on the double porosity mechanism (e.g. Borja & Kavazanjian, 1985; Liu & Borja, 2022), and the isotach model (Bjerrum, 1967; Den Haan, 1994). The latter is based on the isotach rheological principle, which states that the strain rate of the subsurface is uniquely determined by the applied stress. According to this principle, the strain state, defined by the amount of previous compaction and the in situ stress, corresponds to a specific compaction rate that progressively decreases with time. Lines of identical compaction rates can be represented by so-called ‘isotachs’. The model uses three main compaction parameters: the recompression ratio (RR), virgin compression ratio (CR), and viscous compression parameter (C_α), the latter to describe viscous compression (Figure 1). CR characterises elastoplastic deformation under first-time (‘virgin’) loading and is determined in a semi-logarithmic plot from the slope of effective stress versus linear strain, while RR characterises fully elastic deformation and is described by the slope of the unloading curve (Visschedijk, 2010). The C_α parameter describes the rate of compression over time that continues when the loading remains constant, representing the viscous compression behaviour of the subsurface layers. An equal amount of strain happens every log unit of time between two isotachs (Figure 1). Under certain conditions, the viscous compression rate can increase and deviate from this logarithmic path of strain over time, commonly referred to as tertiary compression (e.g. Den Haan & Edil, 1994). This phenomenon poses challenges for modelling, but a consistent inclusion in the isotach model has not been formulated, and hence goes beyond the scope of this article. The isotach model is often used in practise in countries such as the Netherlands, and has been applied successfully in large-scale regional studies, predicting land subsidence in deltaic areas where extensive groundwater extraction leads to significant subsurface compaction (e.g. Bakr, 2015; Keogh et al., 2021; Minderhoud et al., 2017) or where extensive surface drainage leads to compaction of clay and peat in the subsurface (Fokker et al., 2019; Koster et al., 2018; van Asselen, 2010). On a local scale, it has proven valuable in construction projects, such as the design of foundations for water retaining structures and infrastructure, where long-term settlement is a design quality requirement and must be accurately predicted to minimise structural damage by residual settlement (e.g. Kurihara et al., 1994; Nishimoto & Hayashi, 2008; Yamazoe et al., 2025).

Figure 1. Stress-Strain diagram displaying standard compaction behaviour, indicated with the blue line. Elastic compaction occurs along the RR line up to the preconsolidation stress. Beyond the preconsolidation stress compaction occurs along the 1 day isotach. When the final stress has been reached, compaction follows a vertical path in the diagram crossing mutliple isotachs. The rate can be determined with the isotachs, indicated by the red lines. Between every isotach is a factor 10 in time. Isotachs in this diagram are linear, parallel and equidistant. RR: recompression ratio.

Despite its usefulness in application, the isotach model still faces challenges in calibration and validation, particularly in determining C_α. The C_α value of organic subsurface layers shows large variability in laboratory measurements (Den Haan & Kruse, 2007; Fox & Edil, 1996; Mesri & Ajlouni, 2007). When using the isotach model to predict viscous compression, standard values and ranges for C_α are applied that are based on major lithology classes, like peat and clay (NNI, 2016), and on compactness (e.g. Mesri & Castro, 1987; Mesri & Godlewski, 1977). C_α has also been predicted by using a derivative of the pore pressure (Cosenza & Korošac, 2014), further highlighting that choosing an appropriate value for C_α can be difficult. However, uncertainties in the chosen value for the C_α parameter lead to considerable uncertain or incorrect estimates of long-term subsidence. In addition, this approach disregards the heterogeneity of the material observed within these major lithological types, which could be key to understand the large variability seen in organic layers.

Therefore, this study aimed to determine: (1) how viscous compression parameter values vary with material characteristics of peat and clay, and (2) to what extent knowledge of these material characteristics improves the estimation of the viscous compression parameter for different subsurface materials in the isotach model.

To achieve these aims, we analysed an extensive dataset comprising material characteristics and geotechnical parameters, including the viscous compression parameter, from peat, clay and peaty clay/clayey peat samples from the Holocene subsurface of the Rhine-Meuse delta and coastal zone of the Netherlands. The extensive dataset enabled a robust analytical approach, capturing a wide range of lithologies and depositional settings of the material. The full dataset was used to identify lithological classes based on material characteristics, i.e. water content, void ratio, saturated density, dry density and particle density, followed by regression analysis to develop predictive models for C_α based on material characteristics, and to determine to what extent C_α varies within and between lithology type clusters. Together, these steps enabled a detailed evaluation of how variation in material composition influences viscous compression behaviour, supported by the statistical strength and representativeness of the underlying dataset.

Materials and methods

For this study we used a dataset composed by Deltares Research Institute (an applied research institute that works on innovative solutions in the field of water and subsurface) that contains the results of geotechnical laboratory tests of 603 undisturbed peat and clay samples. The samples were collected for research and local advisory projects over multiple years at various locations in the Rhine-Meuse delta, the coastal peatlands of the western Netherlands and in the marine clay area in the Northern Netherlands (Figure 2). Different sampling methods were used, including manual, piston and mechanical coring (pulse drilling, continuous drilling) of which the latter involved also heavy machinery to collect samples with a large diameter of 40 cm and thickness of 80 cm (Deltares large diameter sampler [DLDS]; Zwanenburg, 2017). The dataset entries comprise geotechnical compression test and material characterisation test results. The main texture class of 232 samples in the dataset was labelled as clay and of 349 samples as peat, based on visual description in the field or laboratory during the projects using the NEN5104 lithological texture classification (NNI, 1989). Of the 22 samples the main texture class was not labelled. In the dataset coordinates and absolute or relative depth of the samples were unfortunately predominantly given as averages grouped by sampling project, without giving information on relative positions of the samples, rendering them unusable for further spatial analyses.

Figure 2. Paleogeographic map of the Netherlands for 2000 CE (after Vos et al., 2018). The Pleistocene deposits are depicted in (light) orange colours. Light green areas indicate the diked Holocene floodplains and saltmarshes, dark green areas reclaimed land, and brown areas peatlands. Urban areas are shown in red.

The dataset provides for each sample the results of both incremental loading oedometer tests and constant rate of strain (CRS) tests. Oedometer tests measure the height of a sample over time when subjected to stepwise loading, whereas in CRS tests a constant strain rate is enforced while measuring the load needed to accomplish that as well as the height of the sample over time. Testing procedures for both the oedometer and CRS tests followed from the standard test method for one-dimensional consolidation properties of saturated cohesive soils as described in ASTM D4186-06 and NEN-EN-ISO 17892-5 (ASTM, 2006; ISO, 2017). These procedures prescribe that C_α is determined from the slope of the linear relationship between strain and the logarithm of time during loading steps for which the applied vertical stress exceeds the preconsolidation stress, as is described by Equation 1:

The slope was calculated from strain measurements obtained from 1 day after application of the respective load increment. Consequently, C_α was evaluated only under normally consolidated conditions (overconsolidation ratio = 1). These conditions correspond to a reference time of 1 day and the reference isotach in a stress-strain diagram (Figure 1), which facilitates comparison of C_α values between samples with different loading histories. No noticeable differences were observed between the parameter values of similar samples of the two types of tests. Analysis of the applied stresses, elapsed time, and raw and normalised sample height data produces various compression parameters used in modelling, which were also recorded in the dataset. This includes C_α as the main parameter used in the relation between material characteristics and viscous compression in this research.

Next to measured data, five material characteristics from the dataset were considered in this study: (1) water content by mass (WC), (2) wet bulk density (ρ_wet), (3) dry bulk density (ρ_dry), (4) particle density (ρ_solids), and (5) void ratio (e). These material characteristics were chosen as they describe the solid and fluid components of the subsurface material. Prior to analysis, the dataset was filtered to remove missing or extreme values based on the criteria outlined in Table 1, to ensure the robustness of the clustering results. The cut off value of C_α was set at the maximum value deemed realistic for this parameter; samples with higher C_α values were rejected for analysis. As upper acceptable limit for the particle density reported in the dataset we determined the maximum particle density reported for clay of 2860 kg/m³ (Schjønning et al., 2017). After filtering, a total of 553 samples was used.

Table 1. Selection criteria for the dataset to exclude zero, negative and missing values, and extreme values that may affect the clustering.

Parameter [unit]	Criteria
Viscous compression parameter [–]	0 < Cα < 0.1
Water content [%]	WC > 0
Wet bulk density [kg/m3]	ρwet > 0
Dry bulk density [kg/m3]	ρdry > 0
Particle density [kg/m3]	0 < ρsolids < 2860
Void ratio [–]	e > 0

Firstly, the value range of each material characteristic per main texture class indicated in the dataset was determined. The material characteristics were then checked for collinearity using Pearson’s correlation before applying Principal Component Analysis (PCA) to enable low-dimension clustering and regression. Clustering of the samples based on the PCA results was done using a Gaussian Mixture Model (GMM). The following clusters, or lithological classes as they are called henceforth, show the effectiveness of using these cluster to determine the viscous compression compared to exisiting parameter ranges. Subsequently, multidimensional regression of C_α on a minimal set of material characteristics was applied to both the full dataset as well as to the identified clusters individually. An overview of these steps is presented in Figure 3.

Figure 3. Data processing system diagram showing the various steps taken and the associated data types. Data output from compression tests was split into geotechnical parameters, material characteristics, main texture class and measured data, after which geotechnical parameters and material characteristics were filtered for NA (missing values), zero and unrealistic values. The characteristics were checked for correlation and collinearity and clustering used the output of the PCA. Regression analysis was performed on both the full dataset as well as the clusters with the appropriate regression model. PCA: Principal Component Analysis.

Correlation analysis

Since several material characteristics given in the dataset are mutually related or can be derived from each other (e.g. sample water content from dry and wet bulk densities), to reduce the number of mutually correlated and redundant variables in the clustering and regressions, first a Pearson correlation matrix was constructed to determine the internal linear relationships among the material characteristics. The viscous compression index C_α was only included to assess its correlation with the material characteristics, but not to take into consideration for clustering. Following this, PCA was applied to the material characteristics to capture the greatest variability in the dataset (explanation provided in Appendix 1). Using the principal components instead of the original material characteristics benefited the clustering by separating groups of data points as much as possible. Furthermore, any collinearity between the material characteristics was eliminated with the method.

Organic fraction determination

Naturally, mixtures of clastic and organic material occur in low-lying coastal and delta systems, represented by a lower organic fraction than peat. However, the data cannot directly be separated on organic fraction based on the available material characteristics. To determine the lithological class of the clusters a quantification of the organic and clastic fractions was needed. Since sample-specific sand, silt and clay fractions are unknown, a general clastic fraction was used for the clastic mass density (ρ_clastic). Clay, silt and sand generally have particle mass densities of 2700 kg/m³, 2700 kg/m³ and 2650 kg/m³, respectively, while organic mass density (ρ_org) is 1470 kg/m³ (Erkens, 2009; Poelman, 1975; TAW, 1996). Organic fraction f_org was then calculated using (adjusted from Erkens, 2009; Poelman, 1975; TAW, 1996):

Particle densities (ρ_s) reported in the dataset range between 400 and 3000 kg/m³, of which the extremes are thus lower and higher than the clastic and organic mass densities that make up these samples. Therefore, sample particle densities have been recalculated using the indicated final water content, as full saturation can then be assumed, and final dry bulk density, resulting in a range between 1000 and 2500 kg/m³, apart from a small number of outliers (Appendix 2). These recalculated lower and upper particle density limits were adopted as effective values for ρ_org and ρ_clastic in the organic fraction calculation. This approach assumes that the dataset contains samples approaching end-member compositions, i.e. nearly ‘pure’ peat at the low-density end and nearly ‘pure’ clay at the high-density end. Organic fraction values exceeding 100% were capped at 100%, while negative values were set to 0% to account for all outliers. The organic fraction classification was then defined using the calculated organic fraction and the threshold indicated by the red line in Figure 4, which represents silty clay in terms of clay/silt/sand mixture as the clastic end member. Silty clay was selected as a representative clastic material because samples with higher organic content are predominantly clay-rich, whereas samples with lower organic content exhibit a wider range of clastic compositions.

Figure 4. Classification of samples by their weight percentage of organic material (after De Bakker & Schelling, 1966; NNI, 1989). The red line indicates samples with a clastic component consisting of 42.5% clay, 42.5% silt and 15% sand, which is the composion of silty clay in the USDA classification system (Schoeneberger et al., 2012). The percentages in red denote the organic fraction at the lithological type boundaries.

Cluster analysis

To investigate which groups of samples could be defined by their material characteristics, how these groups link to lithological types, and how viscous compression behaviour varies between different lithological types, a multi-step clustering approach was applied (Appendix 3). Firstly, K-means clustering was applied to the PCA-transformed data, using 25 random initialisations to ensure robustness of the clustering results (Hartigan & Wong, 1979; Appendix 3.1). The output of K-Means served as initialisation for a GMM, enabling the identification of probabilistic cluster memberships (Dempster et al., 1977; Appendix 3.2). As an alternative method, Density-Based Spatial Clustering of Applications with Noise (DBSCAN) was considered (Ester et al., 1996, Appendix 3.3). However, when using the commonly recommended elbow point method and the dimensionality rule to determine parameters, DBSCAN produced numerous small clusters and classified a significant portion of the data as outliers. To integrate the strength of outlier detection observed in DBSCAN, a threshold of 0.99 for cluster probability was applied within the GMM. This threshold excluded datapoints with a lower probability of cluster assignment and instead labelled them as unclassified. Both K-Means and GMM were applied across a range of 2–10 unique clusters. To choose the optimal number of clusters three statistical metrics were obtained to assess the quality of the clustering: the silhouette score, the BIC score and the Jensen-Shannon Divergence (JSD) score (Appendix 4).

Regression analysis

To predict C_α, (multiple) linear regression of C_α was applied on the material characteristics using the entire dataset. In addition, separate regressions were performed for each cluster defined by the earlier clustering analysis. For the full dataset regression, linear, quadratic, logarithmic and non-linear models were tested for the independent variables. Model performance was evaluated based on mean squared error (MSE) and the statistical significance of the regression terms (p-values), while models with a significant amount of negative predicted C_α values, which go against the nature of viscous deformation, were rejected. The regression model was left unchanged for the cluster-specific regression to evaluate the usefulness of the full dataset regression per cluster. This approach allowed for assessing how well the full dataset regression generalised to individual clusters, ensuring model consistency while avoiding unnecessary complexity.

Results

Boxplots of the key material characteristics and viscous compression parameter (C_α) values for clay and peat as indicated in the dataset reveal distinct characteristic differences between the clay and peat types indicated in the dataset (Figure 5). Differences between clay and peat are evident especially in terms of water content and void ratio, where peat exhibits significantly higher values and wider interquartile ranges (IQRs) than clay. Wet bulk densities and dry bulk densities are generally lower for peat with low variability, consistent with its high organic matter content, which has a lower specific density. Finally, particle density is substantially lower in peat, with a skewed distribution and both high-end and low-end statistical outliers, whereas clay shows mainly low-end statistical outliers. Considerable differences between clay and peat are found in the viscous compression behaviour. Clay has a median viscous compression parameter value of 0.0088, with first and third quartiles at 0.0039 and 0.0149, and a maximum value of 0.0682. Peat shows significantly higher values than clay for C_α, with a median of 0.0313 and first and third quartiles of 0.0261 and 0.0373, respectively. Maximum and minimum C_α values for peat are 0.0588 and 0.0043, respectively. The lithological distributions of C_α underline the fundamental differences in physical composition and characteristics between clay and peat. These differences in material characteristics also explain the differences in the viscous compression coefficient values. The higher values of C_α seen for peat can be related to the high water content and void ratios, representing the large pore volume, which can be partially, yet more rapidly compressed through viscous compression. Furthermore, low bulk and particle densities reflect the hollow fibres, which in turn can be compressed, unlike clay particles.

Figure 5. Boxplots showing the range of viscous compression parameter C_α and five material characteristics per main lithological class indicated in the dataset. Dots indicate statistical outliers within the shown class.

Material characteristics correlation

Strong – and obvious – positive correlations were found between water content and void ratio (p = 0.93) (Figure 6a), wet and dry bulk density (p = 0.99), and bulk density and particle density (p > 0.8). The two groups of characteristics (water content/void ratio versus bulk/particle density) showed strong negative correlations with each other (p < -0.7), apart from void ratio and particle density, confirming the expected collinearity between material characteristics. The viscous compression parameter (C_α) positively correlates with water content and void ratio, and strongly negatively to all density characteristics.

Figure 6. (a) Correlogram of five material characteristics and the viscous compression parameter C_α showing internal correlation. The correlation coefficients are ordered using the First Principal Component (FPC) method. Blue colours indicate a positive correlation, while red colours indicate a negative correlation. (b) Principal component analysis result of material characteristics of peat and clay samples. The first two principal components are shown on the x- and y-axes, together explaining 95% of the total variance of the dataset.

The PCA results further reflect these dependencies. The PC1 axis clearly reflects the contrast between the two characteristics groups: high PC1 values are associated with high water content and void ratio, while negative PC1 values were dominated by the high density values. PC1 alone explains 85% of the dataset’s variance, with the first two PCs explaining 95% of total variance (Figure 6b; Appendix 1). All variables have small negative values on PC2, indicating limited orthogonal variance. PC3 explains an additional 4.3% of the variance, accumulating to 99.3% in total.

Given this high degree of collinearity of material characteristics dimensionality reduction was necessary to avoid redundancy and overfitting in clustering. Since principal components eliminate collinearity while retaining nearly all variance, the first three principal components (PC1–PC3) were used for clustering. The rotation matrix used to derive the PC values is included in Appendix 1. Water content and dry bulk density were subsequently selected as representative predictive variables for the regression models, as they (1) capture the core contrast between the two correlated variable groups, (2) are commonly measured in geotechnical investigations in a single, easy to execute test, and (3) showed clear and opposite relationships with C_α.

Sample clusters by material characteristics

After determining nine different configurations of clusters with GMM, the statistical assessments were not unanimous on an optimal number of clusters (Appendix 4). Nonetheless, based on these assessments, four clusters were chosen as most optimal to use for further analysis. In total, 50 samples were not assigned to a cluster as a result of the classification confidence threshold of 0.99. The four clusters show distinct differences in their material characteristics, organic fractions and viscous compression parameter (C_α), (Table 2, Figure 7).

• Cluster 1 is characterised by relatively high bulk and particle densities, along with low water contents and void ratios. It exhibits the lowest organic fractions and the smallest C_α values of all clusters. The distribution of C_α has a particularly narrow IQR and minimal spread between extremes, indicating a consistently low viscous compressibility within this group.
• Cluster 2 has lower bulk and particle densities than Cluster 1, with moderate water contents and void ratios. The organic fraction is higher, and the median C_α is over four times greater than that of Cluster 1, suggesting a clear increase in viscous compression susceptibility. While the IQR of C_α remains moderate, the overall spread is wider, indicating greater internal variation than in Cluster 1.
• Cluster 3 is marked by the lowest particle densities and contains the highest organic fractions. Water content and void ratio values are also high. This cluster exhibits the highest median C_α, with values more tightly distributed than in Cluster 2, though the values are overall higher. These features suggest a strong association between high organic content and increased viscous compressibility.
• Cluster 4 has slightly higher water content and similar wet bulk density as Cluster 3, but has substantially higher void ratios and less variation in particle densities. Although the organic fraction is similarly high, the median C_α is notably lower than in Cluster 3, with least variation of all clusters. This suggests a relatively uniform material behaviour with moderate viscous compression.

Table 2. Statistical summary (median, IQR, and range) of the material characteristics and corresponding lithological classes for each cluster.

		Wet bulk density (kg/m3)	Dry bulk density (kg/m3)	Particle density (kg/m3)	Water content (%)	Void ratio (–)	Organic fraction (%)	Viscous compression parameter (–)	Lithological class
Cluster 1	Median	1627	1056	2568	54	1.347	24.1	0.0074	Slightly Humic Clay – Peaty Clay
	IQR	1523–1769	916–1272	2376–2651	39–68	0.995–1.660	16.7–32.2	0.0035–0.0123
	Min–Max	1291–2308	563–2062	1000–2854	9–140	0.185–2.705	1.6–100.0	0.0002–0.0296
Cluster 2	Median	1163	388	1972	202	3.915	56.0	0.0296	Clayey Peat – Peat
	IQR	1114–1243	313–460	1612–2197	165–251	3.186–4.789	49.3–65.6	0.0224–0.0374
	Min–Max	998–1364	232–910	1075–2677	10–361	0.647–7.044	25.3–100.0	0.0108–0.0682
Cluster 3	Median	990	121	937	727	6.510	85.2	0.0384	Fibrous Peat
	IQR	964–1010	96–150	743–1079	577–898	5.609–7.280	74.2–92.1	0.0338–0.0439
	Min–Max	878–1214	69–260	414–2631	366–1236	3.808–20.463	1.7–100.0	0.0231–0.0588
Cluster 4	Median	968	97	1482	902	14.338	76.7	0.0293	Decomposed Peat
	IQR	950–987	85–115	1453–1512	754–1015	12.075–16.117	68.1–82.0	0.0253–0.0337
	Min–Max	837–1111	61–187	1413–1828	426–1470	7.586–24.620	24.6–100.0	0.0178–0.0584
IQR: interquartile range.

 

Figure 7. Boxplots showing the ranges of five material characteristics, organic fraction and the viscous compression parameter for four clusters assigned with a Gaussian Mixture Model (GMM). The material characteristics serve as input variables for this GMM through their principal components. Outliers outside of any cluster are categorised as NA (missing values) and not included these plots. Organic fraction is added to determine the lithological class of the clusters. Clusters: 1 = Clay, 2 = Clayey Peat, 3 = Fibrous Peat, 4 = Decomposed Peat.

The calculated organic fraction values, in combination with Figure 4, were used to infer the dominant lithological composition of each cluster. As the full range (minimum to maximum) of organic fractions across each cluster spans nearly the entire spectrum of materials shown in Figure 2, the IQR was used as a more representative basis for classification. Based on the IQR, Cluster 1 corresponds to materials classified as clay, slightly to strongly humic clays or peaty clays in Figure 2. Cluster 2 aligns with clayey peat, and peats containing a relatively high proportion of clastic material. Clusters 3 and 4 both represent the highly organic peat class. Compared to Cluster 3, Cluster 4 is characterised by higher particle density, a doubled void ratio, and a modest increase in water content. The slightly higher water content in Cluster 4 is consistent with its higher particle density: if we take hypothetical fibrous and decomposed peat samples with equal dry bulk densities of 100 kg/m³ and particle densities of 950 and 1500 kg/m³, their resulting porosities would be 0.895 and 0.933, respectively. Assuming a water density of 1000 kg/m³, this corresponds to water masses of 895 kg/m³ and 933 kg/m³ per cubic meter of wet soil, giving similar water contents of approximately 895 and 933% (water mass / solid mass). The void ratio, however, shows a much larger difference of 8.52 versus 14.00, reflecting the greater total pore volume in decomposed peat with similar bulk densities. These characteristics are consistent with more humified material, which is typically associated with higher bulk densities and reduced compressibility (Berry, 1983; Huat et al., 2009; Mesri & Ajlouni, 2007; O’Kelly & Pichan, 2013). The differences in median C_α values between Clusters 3 and 4 align with this interpretation as decomposed organic matter is less prone to reorientation and viscous compression due to smaller size of the particles and lower amount of water inside the structure. Given the overall characteristics of the clusters, we interpret and refer to them as the following lithological classes: Cluster 1 will be referred to as ‘Clay’, Cluster 2 as ‘Clayey Peat’, Clusters 3 as ‘Fibrous Peat’ and Cluster 4 as ‘Decomposed Peat’.

Regression models

The dataset shows large variation of the viscous compression parameter, and also shows a considerable range within each lithological class. This introduces uncertainty in determining accurate ‘default’ parameter values for practical modelling cases where no laboratory tests are available. To address this, multiple regression models were developed to predict C_α from material characteristics, both for the full dataset and for the individual lithological clusters identified earlier. Among the material characteristics, water content (WC) and dry bulk density (ρ_d) showed the strongest correlations with C_α under logarithmic transformation, and thus for the full dataset the best-performing regression model was:

The model achieved an R² = 0.61 (p < 0.001), indicating a moderately strong fit. As expected, C_α increases with increasing water content and decreases with increasing dry bulk density (Figure 8a & b). The residual standard error (RSS) was 0.0088 (Figure 9a). Larger residuals were observed for high original C_α values, especially at high water contents and low bulk densities, where predicted values tend to flatten around 0.040. In contrast, original values show larger variability between approximately 0.020 and 0.060. The model tends to overpredict low C_α values. These patterns highlight the difficulty of accurately capturing the viscous compression behaviour with a single regression model, especially of highly organic, water-rich peat layers.

Figure 8. Scatter plots of viscous compression parameter against (a) dry bulk density, entire dataset, (b) water content, entire dataset, (c) dry bulk density per cluster, (d) water content per cluster. The colours for (a) and (b) indicate original and predicted values. The colours for (c) and (d) indicate which cluster a data point belongs to. Circles and diamonds represent the original and predicted values, respectively.

 

Figure 9. Original Cα values vs fitted Cα values. (a) Fitted values based on multiple linear regression using all original data points (R² = 0.61), (b) fitted values based on multiple linear regression per cluster (R² = 0.51, R² = 0.14, R² = 0.05, R² = 0.02 in order of the legend). Indicates viscous compression after a log-unit of time. The total distance covered by the blue and red arrows is schematic representation of the total compaction. The arrows are not to scale.

Regression models were also applied to the individual lithological clusters to assess whether predictions could be improved relative to the full dataset model and to better understand the applicability of the global regression across different material types. All cluster-specific regressions used a logarithmic transformation of both water content and dry bulk density, consistent with the approach for the full dataset. Negative predicted C_α values were set to zero as negative values are physically impossible. These models did not outperform the full regression in terms of predictive accuracy and thus cannot be used as an alternative to the full regression. However, they provided insight into the dominant controls on viscous compression within each cluster (Figure 8c & d). The resulting regression equations are:

Differences between clusters primarily concerned the relative contribution of each predictor. In Clay and Clayey Peat, dry bulk density was the dominant factor influencing C_α, while water content had limited explanatory power especially in Clay. In the Fibrous Peat and Decomposed Peat clusters, both predictors contributed to the model, with dry bulk density generally having a stronger influence.

However, the predictive performance of the cluster-specific models was generally lower than that of the full dataset model. Only the Clay regression showed a moderately good fit, with R² = 0.52 (p < 0.001) and an RSS of 0.0041 (Figure 9b), outperforming the full model within that class – though it occasionally predicted negative C_α values, which are physically unrealistic. The Clayey Peat model performed worse than the full dataset model, with R² = 0.14 (p = 0.004) and RSS of 0.0103. Both Fibrous Peat and Decomposed Peat regressions had very low predictive power (R² = 0.02, p = 0.50 and R² = 0.05, p = 0.01, respectively), and predicted C_α values clustered around 0.040 and 0.030, despite large variation in measured values. These results further support the use of the full dataset regression for general prediction purposes, especially when cluster identity is unknown or when local calibration is impractical.

Alternative regression approaches were tested for the clusters to improve fit for the Clayey Peat and Peat clusters, but none substantially increased predictive performance (R² ≤ 0.20). In summary, the full dataset regression provides a robust and practical tool for predicting the viscous compression parameter across all lithologies. Lithology-specific regression can be used to evaluate the application of the full dataset regresion, rather than provide an improvement in predictive capacity: For Fibrous and Decomposed Peat, full dataset regression values can function as average values, which in combination with ranges may be best suitable for use in modelling and validation.

Discussion

Material characteristics versus C_α

A key point of discussion is whether material classes resulting from this dataset can be used to improve predictions of viscous compression through straightforward preliminary parameterisation, compared to the traditional soil texture classes. Our results separate the dataset samples into more classes than the broad clay–peat division commonly assigned in the field or laboratory, but fewer than the detailed organic soil classes presented in Figure 4. Mainly, the humic subdivision of clay is not represented significantly in their typical material characteristics, which is why the algorithmic clustering approach is unable to distinguish these classes. This underlined by the cluster robustness (Appendix 5), which shows that when five clusters are determined Clay is unchanged and rather Clayey Peat is divided in two. Within peat, however, the clustering provides a particularly meaningful distinction by separating fibrous and decomposed peat. This shows that the state of the organic matter, when it dominates the material, is a critical determinant of viscous compression behaviour.

The distinction between fibrous and decomposed peat is noteworthy because it reveals contrasting compressibility behaviour. While fibrous peat shows higher C_α values, decomposed peat exhibits lower C_α values despite its higher void ratios. This contradicts the common assumption that greater void space leads to greater potential for compaction and viscous compression. Interestingly, C_α values for decomposed peat are more comparable to those of Clayey Peat, even though the two classes differ in most material characteristics. The main consistencies are found in particle density, and to a lesser extent in organic fraction, both of which show patterns aligned with C_α. One plausible explanation is that particle density reflects the mobility of soil particles: denser organic particles, which retain less water within their microstructure, are less prone to rearrangement and therefore less susceptible to viscous compression (Le et al., 2012; van Elderen et al., 2025). This may be because denser particles require a greater initial force to overcome inter-particle friction.

To evaluate the practical relevance of the observed C_α values in our dataset, it is essential to place them in the context of existing reference values used in land subsidence assessments. A comparison with values from Dutch geomechanical practice, listed in NEN-EN 1997-1 (NNI, 2016), and from delta-scale land subsidence modelling in the Mississippi Delta (Keogh et al., 2021, derived from NNI, 2006) provides insight into how our C_α values influence subsidence calculations (Table 3). Although the NEN-standard values were designed to provide lower and upper bounds for settlement calculations rather than to reflect natural variability, they still offer a useful benchmark against which to situate our results. In contrast, our dataset is intended to capture the statistical distribution of C_α values and to explore explanatory trends, thereby offering a broader and more representative frame of reference. Bulk density values correspond relatively well between the NEN-standard and our classes, but the agreement in C_α values is far less consistent. The NEN-standard C_α values for clay, sandy clay, and organic clay classes fall within the range of our Clay class, while the Mississippi Delta C_α values for clay and silty clay loam textures are slightly lower, which might reflect a higher silt percentage compared to the Dutch samples. In contrast, the C_α ranges of Clayey Peat, fibrous peat, and decomposed peat extend well beyond the NEN-standard C_α values, in some cases by a factor three. These differences suggest that our classes provide a more refined framework for identifying compressibility variations within peat layers than is possible using current standard categories alone. They also imply that current assessments may underestimate viscous compression in organic-rich subsurface layers, potentially leading to substantial underestimation of long-term land subsidence.

Table 3. Characteristic values of the viscous compression parameter and wet bulk densities from clustering, and the standard used for geotechnical risk assessments in industry (NNI, 2016) and Mississippi Delta subsidence modelling (Keogh et al., 2021).

Subsurface texture	Clustering		Dutch standard		Mississippi delta
	ρwet (kg/m3)	Cα (–)	ρwet (kg/m3) ± 5%	Cα (–) ± 25%	ρwet (kg/m3)	Cα (–)
Clay	1523–1769	0.0035–0.0123	1400–2000	0.003–0.013	1800	0.0025–0.0046
Sandy clay	-	-	1500–2100	0.0007–0.009	-	-
Silty clay loam	-	-	-	-	1900	0.0020–0.0037
Organic clay	1523–1769	0.0035–0.0123	1300–1600	0.008–0.015	-	-
Clayey peat	1114–1243	0.0224–0.0374	-	-	-	-
Peat (fibrous & decomposed)	950–1010	0.0253–0.0439	1000–1300	0.012–0.023	1050	0.0153–0.0230
Values for the clustering represent the inter quartile range of the lithological class. Values for both the standard and the Mississippi Delta represent the low and high averages. The NNI standard provides an indication of consistency (loose, average, compact) per subsurface texture, which has been merged in this table for the benefit of the comparison with the clustering results. NNI: Nederlands Normalisatie Instituut.

A key question left unanswered is how much of the observed compaction originates from viscous compression, as this determines the broader interpretation of the results. One way to assess the relative contribution of viscous versus virgin compression is by examining their ratio in the isotach model. Originally, Mesri and Godlewski (1977), and later extended by Mesri and Castro (1987), defined this ratio as C_c/C_αe, where C_c is the coefficient of compression and C_αe the coefficient of secondary compression, both based on the void ratio – logarithm of stress curve. These parameters correspond to CR and C_α, respectively, which are based on the strain – logarithm of stress curve, and the relationship remains valid (Den Haan & Kruse, 2007). A high CR/C_α ratio indicates that virgin compression accounts for a larger share of total compressibility, while lower ratios reflect a stronger influence of viscous compression. Importantly, the ratio represents relative susceptibility rather than absolute dominance of one process over the other.

To illustrate this, we take the mean values for Clay and Clayey Peat: if C_α increases from 0.008 to 0.031 while CR increases from 0.18 to 0.40, the ratio falls from 23 to 17. Although viscous compression contributes relatively more in the second case, virgin compression still accounts for a substantial absolute share of total compaction. The contribution in decomposed and fibrous peat is derived in the same way and results in ratios of 17 and 13. Typical ratios reported by Mesri and Godlewski (1977) are ~50 for granular soils, ~25 for inorganic clays, ~20 for organic clays, and ~17 for peat. Den Haan and Kruse (2007) further proposed an S-shaped curve describing how CR/C_α evolves with increasing wet bulk density (Figure 10). Applying this ratio to our dataset reveals two main observations (Figure 10). Firstly, samples within the Clay cluster tend to show significantly higher CR/C_α ratios than predicted by the reference curve, implying that these materials are less prone to viscous compression than according to Den Haan and Kruse (2007). Secondly, decomposed peat diverges from the expected pattern of decreasing ratios with decreasing wet bulk density. Moreover, while decomposed peat and clayey peat show similar C_α values, a comparable relationship is observed between fibrous peat and clayey peat for the CR/C_α ratio. The higher ratios of decomposed peat align with the lower mean C_α values for this class compared to fibrous peat at similar bulk densities. Two hypotheses may explain this behaviour. One is that decomposed peat particles, being denser and containing less internal water, are less prone to rearrangement, so more of the compaction occurs through virgin compression of macropores. An alternative explanation is that once a material is more susceptible to virgin compression, this limits the extent to which viscous compression can develop. Testing this hypothesis would require consolidation experiments under increasing stresses with intermittent phases of viscous compression.

Figure 10. Ratio between virgin compression parameter CR and viscous compression parameter Cα as a function of dry bulk density. The black curved line indicates the relationship found by Den Haan and Kruse (2007) on samples from Sliedrecht, the Netherlands. CR: compression ratio.

Practical application

This study demonstrates that regression across the full range of material characteristics can provide meaningful estimates of the viscous compression parameter to use in modelling when laboratory determination is in process or unavailable. However, predictions at the material class level remain subject to uncertainty, particularly for peat. In highly organic samples with low dry bulk density and high water content, variability remains large, suggesting that simple empirical models should be applied cautiously. For practical purposes, class-based representative values are therefore recommended. For clayey peat, fibrous peat, and decomposed peat, upper and lower limits of 0.060 and 0.020 for C_α, respectively, can be used as best- and worst-case scenarios to guide design and assessment. Relative to the median C_α values of each material, taken as representing 100% viscous compression per logarithmic time unit, the upper limit corresponds to approximately 203, 156, and 205%, while the lower limit corresponds to 68, 52, and 68% for clayey peat, fibrous peat, and decomposed peat, respectively.

Future research

This study can serve as starting point for prediction of other geotechnical model parameters, for example, the virgin CR or the RR. With this dataset, one similar to it, or using the material characteristic ranges resulting from this study one can determine mean values for these parameters and look into regression estimation similar to the method presented here. Another improvement can be made by incorporating more data into the analysis, for example, from other countries than the Netherlands, to validate the usage of the results in various conditions and potentially discover new relationships between material characteristics and viscous compression. Furthermore, our results spike the question which changes in viscous compression occur when organic subsurface layers undergo decomposition, as we indicate a clear distinction between the viscous compression of fibrous and decomposed peat. This might have a critical impact on the long-term settlement assessments for mitigation options. If C_α values for peat are higher than currently assumed, but decrease over timescales of decades due to decomposition, the viscous compression occuring over a 1 year period will be larger than previously predicted, up to twice the amount for normally consolidated peat using the values from Table 3. Contribution of viscous compression over a 10 or 100 year period might be equal or lower than presumed in current assessments because of the reduction of C_α with decomposition and, therefore, viscous compression rate over logarithmic time decreases rather than staying constant.

Conclusion

Our research aimed to determine lithology-based variation of viscous compression parameter C_α values used in compaction models for land subsidence assessments, and to evaluate whether simple material characteristics can help the determination of viscous compression parameter values for different soft substrate lithologies. An extensive dataset of compression tests on peat and clay enabled the identification of four distinct lithological classes based on water content, wet and dry bulk densities, particle density, and void ratio: Clay, Clayey Peat, Fibrous Peat, and Decomposed Peat. Median viscous compression parameter (C_α) values increase with increasing organic fraction from 0.0074 for Clay to 0.0296, 0.0384 and 0.0293 for Clayey Peat, Fibrous Peat and Decomposed Peat, respectively. The lithological class-specific C_α median values and ranges are an improvement compared to existing standard values. Prediction of the viscous compression parameter can be done relatively well when using the entire dataset (R² = 0.61, p < 0.001) as opposed to fits per individual lithological class. Predictions for Clay are still significant (R² = 0.52, p < 0.001), whereas predictions for Clayey Peat and Fibrous Peat & Decomposed Peat show low significance (R² < 0.15, p ≥ 0.01) with near horizontal fits. These fits provide an adequate average, which should be used in combination with minimum and maximum C_α values of 0.020 and 0.060, corresponding to approximately halving and doubling of the amount of viscous compression over time. Lastly, we highlight that an adaptable viscous compression parameter for peat with shifting material characteristics over time, e.g. due to decomposition, can advance modelling precision.

Acknowledgements

We are grateful for the input of Hans Middelkoop (Utrecht University) during the internal revision process, as the discussion of material characteristics and their interaction improved the quality of the review. Furthermore, we want to thank Deltares research institute for giving access to their geotechnical parameter dataset, which serves as the spine of this article. The research presented is part of the project Living on soft soils: subsidence and society (grantnr. NWA.1160.18.259). This project is funded by the Dutch Research Council (NWO-NWA-ORC), Utrecht University, Wageningen University, Delft University of Technology, Ministry of Infrastructure & Water Management, Ministry of the Interior & Kingdom Relations, Deltares, Wageningen Environmental Research, TNO Geological Survey of The Netherlands, STOWA, Water Authority: Hoogheemraadschap de Stichtse Rijnlanden, Water Authority Drents Overijsselse Delta, Province of Utrecht, Province of Zuid-Holland, Municipality of Gouda, Platform Soft Soil, Sweco, Tauw BV, NAM.

References

ASTM D4186/D4186-06, 2006: Standard test method for onedimensional consolidation properties of saturated cohesive soils using controlled-strain loading. West Conshohocken, PA: ASTM International.

Augustesen, A., Liingaard, M. & Lade, P.V., 2004: Evaluation of time-dependent behavior of soils. International Journal of Geomechanics 4(3): 137–156. DOI: 10.1061/(ASCE)1532-3641(2004)4:3(137)

Bakr, M., 2015: Influence of groundwater management on land subsidence in deltas: a case study of Jakarta (Indonesia). Water Resources Management 29(5): 1541–1555. DOI: 10.1007/s11269-014-0893-7

Berry, P.L., 1983: Application of consolidation theory for peat to design of a reclamation scheme by preloading. Quarterly Journal of Engineering Geology 16(9): 103–112. DOI: 10.1144/GSL.QJEG.1983.016.02.03

Bjerrum, L., 1967: Engineering geology of Norwegian normally-consolidated marine clays as related to settlements of buildings. Geotechnique 17(2): 83–118. DOI: 10.1680/geot.1967.17.2.83

Borja, R.I. & Kavazanjian, E., 1985: A constitutive model for the stress–strain–time behaviour of ‘wet’ clays. Geotechnique 35(3): 283–298. DOI: 10.1680/geot.1985.35.3.283

Chai, J.C., Jia, R. & Hino, T., 2012: Anisotropic consolidation behavior of Ariake clay from three different CRS tests. Geotechnical Testing Journal 35(6). DOI: 10.1520/GTJ103848

Cosenza, P. & Korošak, D., 2014: Secondary consolidation of clay as an anomalous diffusion process. International Journal for Numerical and Analytical Methods in Geomechanics 38(12): 1231–1246. DOI: 10.1002/nag.2256

de Bakker, H. & Schelling, J., 1966: Systeem van bodemclassificatie voor Nederland. De hogere niveaus. Grondboor & Hamer 20(5): 229–229.

Deltares, 2018: D-settlement embankment design and soil settlement prediction: user manual. Delft: Deltares.

Dempster, A.P., Laird, N.M. & Rubin, D.B., 1977: Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society, Series B 39(1): 1–38. DOI: 10.1111/j.2517-6161.1977.tb01600.x

Den Haan, E.J., 1994: Stress-independent parameters for primary and secondary compression. In: International conference on soil mechanics and foundation engineering. Rotterdam: A.A. Balkema, pp. 65–70.

Den Haan, E.J. & Edil, T.B., 1994: Secondary and tertiary compression of peat. In: International workshop on advances in understanding and modelling the mechanical behaviour of peat. Rotterdam: A.A. Balkema, pp. 49–60.

Den Haan, E.J. & Kruse, G.A.M., 2007: Characterisation and engineering properties of Dutch peats. In: Tan, T.S., Phoon, K.K., Hight, D.W. & Leroueil, S. (Eds.): Characterisation and Engineering Properties of Natural Soils. Vol. 3. Proceedings of the Second International Workshop on Characterisation and Engineering Properties of Natural Soils, Singapore, 29 November–2 December 2006, A.A. Balkema, pp. 2103–2133.

Edil, T.B. & Dhowian, A.W., 1979: Analysis of long-term compression of peats. Geotechnical Engineering 10(2): 159–178.

Erkens, G., 2009: Sediment dynamics in the Rhine catchment: quantification of fluvial response to climate change and human impact. Utrecht: Utrecht University.

Ester, M., Kriegel, H., Sander, J. & Xiaowei, X., 1996: A density-based algorithm for discovering clusters in large spatial databases with noise. In: Proceedings of the Second International Conference on Knowledge Discovery and Data Mining, Portland, Oregon. Institute for Computer Science, University of Munich, pp. 226–231. DOI: 10.5555/3001460.3001507

Fokker, P.A., Gunnink, J.L., Koster, K. & de Lange, G., 2019: Disentangling and parameterizing shallow sources of subsidence: application to a reclaimed coastal area, Flevoland, the Netherlands. Journal of Geophysical Research: Earth Surface 124(5): 1099–1117. DOI: 10.1029/2018JF004975

Fokker, P.A. & Orlic, B., 2006: Semi-analytic modelling of subsidence. Mathematical Geology 38(5): 565–589. DOI: 10.1007/s11004-006-9034-z

Fox, P.J. & Edil, T.B., 1996: Effects of stress and temperature on secondary compression of peat. Canadian Geotechnical Journal 33(3): 405–415. DOI: 10.1139/t96-062

Geertsma, J., 1973: Land subsidence above compacting oil and gas reservoirs. Journal of Petroleum Technology 25(6): 734–744. DOI: 10.2118/3730-PA

Hartigan, J.A. & Wong, M.A., 1979: Algorithm AS 136: a K-means clustering algorithm. Applied Statistics 28: 100–108. DOI: 10.2307/2346830

Huat, B.B., Asadi, A. & Kazemian, S., 2009: Experimental investigation on geomechanical properties of tropical organic soils and peat. American Journal of Engineering and Applied Sciences 2(1): 184–188. DOI: 10.3844/ajeas.2009.184.188

ISO 17892-5, 2017: Geotechnical investigation and testing-laboratory testing of soil-part 5: incremental loading oedometer test. Brussels: European Committee for Standardization.

Keogh, M.E., Törnqvist, T.E., Kolker, A.S., Erkens, G. & Bridgeman, J.G., 2021: Organic matter accretion, shallow subsidence, and river delta sustainability. Journal of Geophysical Research: Earth Surface 126(12): e2021JF006231. DOI: 10.1029/2021JF006231

Koster, K., Stafleu, J., Cohen, K.M., Stouthamer, E., Busschers, F.S. & Middelkoop, H., 2018: Three-dimensional distribution of organic matter in coastal-deltaic peat: implications for subsidence and carbon dioxide emissions by human-induced peat oxidation. Anthropocene 22: 1–9. DOI: 10.1016/j.ancene.2018.03.001

Kurihara, N., Isoda, T., Ohta, H. & Sekiguchi, H., 1994: Settlement performance of the central Hokkaido expressway built on peat. In: Advances in understanding and modelling the mechanical behaviour of peat: Proceedings of the internation workshop, 16–18 June, 1993, Delft: A.A. Balkema, pp. 361–367.

Le, T.M., Fatahi, B. & Khabbaz, H., 2012: Viscous behaviour of soft clay and inducing factors. Geotechnical and Geological Engineering 30: 1069–1083. DOI: 10.1007/s10706-012-9535-0

Liu, Y. & Borja, R.I., 2022: Time scales in the primary and secondary compression of soils. International Journal for Numerical and Analytical Methods in Geomechanics 46(8): 1383–1408. DOI: 10.1002/nag.3350

Mesri, G., 1973: Coefficient of secondary compression. Journal of the Soil Mechanics and Foundations Division 99(1): 123–137. DOI: 10.1061/JSFEAQ.0001840

Mesri, G. & Ajlouni, M., 2007: Engineering properties of fibrous peats. Journal of Geotechnical and Geoenvironmental Engineering 133(7): 850–866. DOI: 10.1061/(ASCE)1090-0241(2007)133:7(850)

Mesri, G. & Castro, A., 1987: Cα/Cc concept and K 0 during secondary compression. Journal of Geotechnical Engineering 113(3): 230–247. DOI: 10.1061/(ASCE)0733-9410(1987)113:3(230)

Mesri, G. & Godlewski, P.M., 1977: Time-and stress-compressibility interrelationship. Journal of the Geotechnical Engineering Division 103(5): 417–430. DOI: 10.1061/AJGEB6.0000421

Minderhoud, P.S., Erkens, G., Pham, V.H., Bui, V.T., Erban, L., Kooi, H. & Stouthamer, E., 2017: Impacts of 25 years of groundwater extraction on subsidence in the Mekong delta, Vietnam. Environmental Research Letters 12(6): 064006. DOI: 10.1088/1748-9326/aa7146

Mitchell, J.K. & Soga, K., 2005: Fundamentals of soil behavior (Vol. 3). New York, NY: John Wiley & Sons, p. 558.

Nederlands Normalisatie Instituut (NNI), 1989: Geotechniek: Classificatie van onverharde grondmonsters, NEN 5104. Delft: NNI, 24 pp.

Nederlands Normalisatie Instituut (NNI), 2006: Geotechnics – TGB 1990 – basic requirements and loads. Delft: NEN.

Nederlands Normalisatie Instituut (NNI), 2016: NEN-EN 1997-1: Eurocode 7, national annex The Netherlands; Geotechnical design of structures – part 1: general rules. Delft: NNI.

Nishimoto, S. & Hayashi, H., 2008: A case study of the long-term settlement and life cycle cost of a road on peaty ground. In: Ellis, E., Yu, H.-S., McDowell, G., Dawson, A.R. & Thom, N. (Eds.): Advances in transportation geotechnics. London: CRC Press, pp. 411–414.

O’Kelly, B.C. & Pichan, S.P., 2013: Effects of decomposition on the compressibility of fibrous peat – a review. Geomechanics and Geoengineering 8(4): 286–296. DOI: 10.1080/17486025.2013.804210

Poelman, J.N.B., 1975: Dichtheid van de vaste delen van rivierkleigronden. Boor en Spade XIX: 32-38. Wageningen: Stichting voor Bodemkartering.

Schjønning, P., McBride, R.A., Keller, T. & Obour, P.B., 2017: Predicting soil particle density from clay and soil organic matter contents. Geoderma 286: 83–87. DOI: 10.1016/j.geoderma.2016.10.020

Schoeneberger, P.J., Wysocki, D.A. & Benham, E.C. (Eds.), 2012: Field book for describing and sampling soils. Government Printing Office. Lincoln, NE: Natural Resources Conservation Service, USDA, National Soil Survey Center.

Technische Adviescommissie voor de Waterkeringen (TAW), 1996: Technisch rapport geotechnische classificatie van veen. Report. Delft: Ministry of Transport and Water Management, 102 pp.

Terzaghi, K., 1941: Undisturbed clay samples and undisturbed clays (No. 16). Boston: Harvard University.

van Asselen, S., 2010: Peat compaction in deltas: implications for Holocene delta evolution. Utrecht: Utrecht University.

Van Asselen, S., Stouthamer, E. & Van Asch, T.W., 2009: Effects of peat compaction on delta evolution: a review on processes, responses, measuring and modeling. Earth-Science Reviews 92(1–2): 35–51. DOI: 10.1016/j.earscirev.2008.11.001

van Elderen, P., Erkens, G., Zwanenburg, C., Middelkoop, H. & Stouthamer, E., 2025: Viscous compression of clay and peat. Earth-Science Reviews 260: 104993. DOI: 10.1016/j.earscirev.2024.104993

Visschedijk, M., 2010: Isotachen berekeningen op een sigarendoosje. Geotechniek 14(3): 30.

Xotta, R., Zoccarato, C., Minderhoud, P.S.J. & Teatini, P., 2022: Modeling the role of compaction in the three-dimensional evolution of depositional environments. Journal of Geophysical Research: Earth Surface 127(9): e2022JF006590. DOI: 10.1029/2022JF006590

Yamazoe, N., Nishimura, S., Tanaka, H., Ogino, T. & Kochi, T., 2025: Long-term settlement behavior of peat after unloading and applicability of isotach law. Soils and Foundations 65(1): 101560. DOI: 10.1016/j.sandf.2024.101560

Zhao, D., Gao, Q.F., Hattab, M., Hicher, P.Y. & Yin, Z.Y., 2020: Microstructural evolution of remolded clay related to creep. Transportation Geotechnics 24: 100367. DOI: 10.1016/j.trgeo.2020.100367

Zwanenburg, C., 2017: The development of a large diameter sampler. In: Proceedings of the 19th International Conference on Soil Mechanics and Geotechnical Engineering, Seoul, 17–21 September 2017.