Experimental observation and constitutive modelling of the shear strength of a natural unsaturated soil

. Increasing frequency and intensity of extreme weather events in the Netherlands is raising attention on the unsaturated response of geo-infrastructures, promoting research projects to provide an overview of the impact of unsaturated conditions on the response of shallow soil layers and embankments, and to better address maintenance and mitigation measures. As part of this effort, we discuss the results of standard laboratory tests performed on initially unsaturated samples retrieved from the field and tested in natural conditions, as well as after controlled drying and wetting. The variation of the “undrained” (i.e. at constant water content) shear strength with the degree of saturation obtained from the laboratory tests aligns well with CPT measurements performed in the field. An elastic-plastic constitutive model with mixed isotropic-rotational hardening developed for saturated soft soils was extended to unsaturated conditions by following a robust approach previously developed for compacted clayey soils. Coupling between the mechanical and the hydraulic behaviour is provided by the water retention curve. The model nicely captures the response observed in the laboratory, until extreme dry conditions, which possibly alter the structure of the soil, the peak stress, and the brittleness after failure. The model is capable of reproducing the effects of the previous hydraulic history on the stress-strain behaviour observed from the laboratory tests over a wide range of degree of saturation.


Introduction
In recent years, the Netherlands has experienced several heat waves, hot summers, and temperatures peaks above average [1], with the heat waves in June and July 2019 exceeding 40°C. These hot summers raised attention on the role of unsaturated soil behaviour on the response of geo-infrastructures. Typically, these are assessed based on the results of cone penetration tests (CPT) performed from spring to autumn. Although over this period the upper soils are likely to be unsaturated, the influence of reduced water content on the interpretation of the test was seldom investigated, which hinders reliable quantification of the strength over time.
In a research project initiated by the General Directorate of the Ministry of Public Works, CPTs were performed every 2 to 4 weeks over an entire year at two sites, where in-situ sensors were installed to monitor water content, suction, and pore water pressure. In parallel, an experimental programme was conducted at TU Delft, to investigate the "undrained", constant water content, unsaturated shear strength. The results obtained in the laboratory are presented and analysed with the aid of an advanced constitutive model for soft soils recently developed and extended to account for unsaturated conditions. * Corresponding author: c.chao@tudelft.nl

Field data
The investigation was carried out from the crest at the Maasdijk near Oijen as displayed in Fig. 1. The dyke, built in the 1950s, is made of clay and sandy clay. The subsoil at the location consists of 1.0 to 3.0 m thick Holocene clay layers and sandy clay layers on top of the Pleistocene sand layer [2]. The position of the sensors installed on site is also displayed in Fig. 1, together with the position of the sampling tubes from borehole B002 tested in the laboratory. The reference for the elevations adopted here is the NAP (Normaal Amsterdam Peil).
Cone penetration tests were repeated over time to quantify the variation of the CPT measurements during different periods, both in the dry and the wet seasons [3]. On each date, two series of CPT were performed. Fig. 2 displays the tip resistance, , from one of the two series over the period ranging from 18/09/2019 to 11/05/2020. The daily precipitation records from the station at Megen (station code 903, KNMI) and at Volkel, (station code 375, KNMI) a few km from Oijen are displayed in Fig. 3. As shown in the figure, September 2019 to December 2019 corresponds to a transition from a very dry condition to a very wet condition. On the contrary, from March 2020 to May 2020 the weather data reflects a transition between wet and average dry conditions.   As shown in Fig. 2, the cone resistance in the first two meters, from 7 to 9 m NAP, presents remarkable variations during the monitoring period. The profiles of the cone resistance in the field show a dramatic reduction (5 MPa to 1 MPa) occurred in the first 2 m of soil after the intense precipitation in October 2019. Contrarily, the drying period experienced in the spring of 2020 contributed to an increase in the cone resistance from 1 MPa to values above 4.5 MPa.
The measurements of the pore water pressure from the tensiometers installed at the Maasdijk near Oijen are reported in Fig. 4. The intense precipitation in late autumn and winter 2019 reduced the suction almost to zero in the upper part of the soil (S1 and S2). The tensiometers S1 and S2 showed a noticeable increase in suction starting in the spring of 2020. The most surficial tensiometer, 1 m below the ground, reached a maximum suction of about 180 kPa in early June, while the sensor S2 at 7.6 m NAP experienced cavitation at about 80 kPa later on in August. The deepest sensors S3 and S4 did not show significant suction development over the entire period.

Laboratory tests
Triaxial tests were conducted on undisturbed samples collected from 7.8 m to 9.1 m NAP (Tube 1, 2, 3, 4) in the very wet period at the end of October 2019. To guarantee a sufficiently homogeneous distribution of water content within the sample during drying and wetting, the triaxial tests were conducted on 38 mm diameter samples with a height-to-diameter ratio equal to 2.2. The representative confining stress was set to 15 kPa. Table 1 lists the tests conducted on the samples. Sample T1S1 was first dried under controlled temperature and relative humidity, at T = 17 °C and RH = 70%, for three days, before mounting in the triaxial cell, while sample T3S1 was dried for 6 days. The evolution of water content during the drying stage is reported in Fig. 5.
The plastic and liquid limits of the test samples are about = 0.22 and = 0.35, with an organic content of 4% and specific gravity ranging from 2.61 to 2.67. The particle size distribution for each sample is reported in Fig. 6, and relevant information on the samples at the start of the tests is reported in Table 2.   All the shear stages started from isotropic stress. The results are presented using the deviatoric stress, q, and the axial strain, . The data have been elaborated assuming a cross-sectional area correction of an equivalent cylinder, and membrane correction has been applied to the radial stress as proposed by [4] and [5]. For samples developing a localised shear band at failure, the stress strain relationship is considered representative of the sample behaviour up to the maximum deviatoric stress. The deviatoric stress strain response of the four samples is compared in Fig. 7. The saturated sample T4S2 showed a ductile asymptotic deviatoric response and a barrel deformation mode. In contrast, the deviatoric response turned to brittle as the degree of saturation decreased. At decreasing water content, failure occurred on localised shear planes (samples T4S1 and T1S1) and eventually axial splitting (sample T3S1).    WP4C® [6][7] starting from fully saturated conditions. A shrinkage test was completed on sample T2S3 to estimate the degree of saturation. The sample, 38 mm in diameter and 38 mm in height, was dried at T =17°C and RH = 63% for several days while measuring its mass and dimensions at regular intervals. The shrinkage curve of sample T2S3 is displayed in Fig. 9 in the void ratio, , water ratio, = , plane. The modified van Genuchten's model from [8], ( ) is a correction function. In the equations above, , m are related to the air-entry value and residual water ratio of the soils. Parameter n controls the slope of the inflection point of the curve. The best fitting parameters in Table 3 were adopted to reproduce the data as shown in Fig. 10.

Numerical simulations
To provide better insight into the unsaturated response of the clay for engineering purposes, an elastic-plastic model with mixed isotropic-rotational hardening, recently developed for soft organic soils at TU Delft, was extended to unsaturated conditions following a robust approach previously introduced and discussed [9]. The constitutive equations developed with reference to saturated conditions are based on the approach by Dafalias and coworkers [10]. The yield surface is taken from [11] to provide more flexibility in its shape, through the parameters and , though keeping the formulation rather simple. Written in terms of triaxial variables, the yield surface reads: where = − , is the deviatoric stress, the total isotropic stress is = ( + 2 ) 3 ⁄ , and the isotropic stress acting on the soil skeleton (which becomes the effective stress when the soil becomes saturated) is defined as ^= + and represents the current rotation of the yield surface in the meridian plane. ^0 is the preconsolidation pressure. The plastic potential is expressed with triaxial variable as: where is the critical state stress ratio and ^ is a dummy variable defining the size of plastic potential.
The evolution of anisotropy in the model is described by the rotational hardening rule adopted by [10]: In equation (4) ⟨ ⟩ is the plastic multiplier, and are the slopes of the isotropic normal compression line and recompression line, ^ (^=⁄ ) is the current stress ratio, and is a model parameter controlling the rate of evolution of and controls the target value of for a given stress ratio. Extension to unsaturated conditions is achieved by introducing ^ as a convenient measure of the "effective" isotropic stress, and by making the hardening rules dependent on the degree of saturation [12]. Coupling between the mechanical and the hydraulic response is provided by the water retention model, where the degree of saturation is adopted as a convenient measure of the soil water content. The evolution of the preconsolidation pressure adopted here reads: where 1 and 2 are model parameters, which measure the sensitivity of the preconsolidation pressure to changes in saturation and define the shape of the so- called "loading-collapse" curve. For a detailed description of the conceptual derivation and calibration of this class of models, reference can be made to [12].
The parameters describing the behaviour of the saturated soil were calibrated on the results of the test performed on Tube 12 (4.56 -4.86 m NAP), retrieved below the lowest position of the phreatic surface (at 5.3 m NAP).
The critical stress ratio was defined based on consolidated undrained triaxial compression tests on Tube 12 which gave a friction angle ′ = 31 . The slope of the isotropic normal compression line, unloading-reloading lines, and Poisson's ratio are chosen from a large data set of silty clay from a parallel study [13].
The coefficient of lateral earth pressure at rest for normally consolidated conditions was calculated from Jaky's simplified relationship which gives 0 = 0.48. The initial stress state ( ′ 0 , 0 ) was defined with reference to an average depth of 5 meters depth below the ground surface and a bulk unit weight of 18 kN/m 3 . The inclination of the yield surface, , was chosen to reproduce normally consolidated saturated conditions, = 0.454. The model parameters used in the simulations are reported in Table 4. The preconsolidation pressure corresponding to the initial stress is derived from equation (2), which gives ′ 0 = 63 kPa. All the parameters were calibrated on the saturated soil behaviour, except 1 and 2 , which were tentatively given a value based on the experience on compacted clays, and eventually refined on the current data. The most crucial aspect of the simulation is the proper initialisation of the stress history. The numerical simulation of one of the tests used for calibration, assuming the sample started from normally consolidated conditions is shown in Fig. 11 with the blue line. The simulation clearly underestimates the maximum deviator stress reached by the sample in the test. The experimental stress path followed by the sample suggests that the higher strength reached in the test may come from the history of overconsolidation experienced during the construction of the dyke. A new simulation was run, where the preconsolidation pressure was increased to ′ 0 = 152 kPa to simulate compaction stress of about 90 kPa during construction. For the sake of simplicity, it was assumed that the compaction occurred under 0 conditions. The new simulation, reported with the red line in the same figure shows much better agreement with the experimental results.
Given the previous results, the equivalent compaction stress was added at the beginning of the simulation of the triaxial tests, which were run after imposing the drying history experienced by the samples in the laboratory. The model predictions are presented together with the experimental results in Fig. 12, the model captures the response observed in the laboratory, until extreme dry conditions, which probably change the structure of the soil leading to a brittle response.

Discussion and conclusions
On average, the model predictions show very good agreement with the experimental data for degrees of saturation above 0.50. For lower degree of saturations, the model continues predicting a ductile response at failure, which is reached on the critical state line. On the contrary, the experimental results show a brittle response of the soil after reaching a clear peak stress, which increases with suction, similar to highly overconsolidated soils.
The reason behind the difference with the prediction can be better understood by looking at the stress path predicted by the simulations in Fig. 13 overpressure is generated and failure is reached in almost "drained" conditions. At lower water contents, the stress path starts following the "drained" stress path at the beginning, however, at some stress level it starts bending to the left, indicating an increase in pore pressure -decrease in suction -which still brings the soil to failure on the critical state line. In no case, among the simulated ones, the stress path is able to overpass the critical state line reaching higher strength that the critical one.
Notwithstanding the qualitative differences, which are currently under study, the theoretical prediction of the strength at constant water content, obtained with the advanced model, is satisfactory for a wide range of suctions. The results suggest that the model can be used to preliminarily infer the variation over time of the shear strength due to seasonal climate fluctuations.
To corroborate this conclusion, Figure 14 reports the comparison between the strength predicted by the model and selected data from the field, derived from CPT data near the water content sensors in the summer and autumn of 2019. The undrained shear strength is estimated from the CPTs using a transformation factor equal to 15 taken from [3]. The theoretical prediction is satisfactory and has the advantage of being almost always on the safe side. At present, experimental tests are run at TU Delft with a new setup fully instrumented for unsaturated conditions. The experimental tests will improve the understanding of the unsaturated behaviour of the Dutch highly organic soft soils, to improve future assessment of geotechnical infrastructure under increasing climatic stresses.
The financial contribution of Deltares under the contract "Overeenkomst van Opdracht Laboratorium testen op kleimonsters voor onderzoek sterkte van grond in de onverzadigde zone" 11205262-010 is gratefully acknowledged. In-situ Laboratory (ultimate) Laboratory (peak) Model