Viscoplastic cyclic degradation model for soft natural soils

Cyclic loading affects the long-term response of geostructures build on natural soils. Soft soils are particularly susceptible to the development of large deformations, induced by the repetitive nature of loading. A new vis-coplastic cyclic accumulation model is presented, which is an hierarchical extension of the Creep-SClay1S model, to model the long-term permanent deformation resulting from undrained cyclic loading of natural soft clays. The cyclic cumulative strains are incorporated by means of an additional viscoplastic multiplier. This cyclic visco-plastic multiplier adds four additional model parameters that are derived from undrained cyclic triaxial tests. The model is calibrated using experimental data from undrained cyclic triaxial tests performed on high quality block samples of natural Ons ø y clay, at different average shear stresses, shear stress amplitudes and loading periods. The accuracy of the proposed model is demonstrated by comparing the element level simulations with the experimental data. The applicability of the proposed model is further illustrated with a boundary value problem, where an embankment submitted to cyclic loading is simulated. The use of the new model enables the simulation of the response of cyclic loaded foundations on soft soils, where the serviceability limit state over a long period of time is governing the design.


Introduction
Cyclic loading is of importance for geotechnical structures, when a significant variation of the load in time is expected, such as is the case for offshore foundations, wind turbines, railway tracks, pipelines and natural slopes (Abdelkrim et al., 2003;Prisco et al., 2003;Andersen, 2009).These loads often can be considered quasi-static, as in most cases the loading period is sufficiently long to disregard the inertial effects in the soil caused by the loading.The non-linear nature of behaviour of soils that lead to irreversible deformations that accumulate over time, however, cannot be as easily disregarded.
The current state of the art for quasi-static analysis of monotonic loads, such as the bearing capacity of foundations, includes complex non-linear models that in addition to the equilibrium are able to predict, quite accurately, the deformations.However, this accuracy is not yet obtained for problems involving cyclic loading.When complex constitutive models, specifically designed to capture each loading cycle [e.g.] ( Mroz et al., 1978;Li and Meissner, 2002;Dafalias and Manzari, 2004;Tafili and Triantafyllidis, 2020), are used to solve cyclic problems with a large number of loading cycles, O(10 5 ), issues related to the computational cost and numerical convergence arise (Niemunis et al., 2005).
In order to overcome these difficulties, several authors have proposed empirical models to predict the cumulative cyclic deformation [e.g.] (Li and Selig, 1996;Chai and Miura, 2002;Ren et al., 2018).In general, the empirical models are based on relationships between the cumulative plastic strain and the number of loading cycles (some empirical models also consider stress and soil type).The main advantage of the empirical models is their ease of use, and the limited experimental data required for the soil to make a prediction.However, this is also their biggest limitation, as often their accuracy is poor, as the soil behaviour captured is overly simplified.
An alternative to overcome the limitations of both complex constitutive and empirical models, are the cyclic accumulation models that have emerged in the past years.The working principle of the accumulation models is the combination of a mechanical constitutive relation to model the first load cycle(s), with accumulation functions to predict the deformation for the remaining number of loading cycles (Suiker and de Borst, 2003;Niemunis et al., 2005).The accumulation models were initially developed for granular materials (Suiker and de Borst, 2003;Niemunis et al., 2005;Karg et al., 2010;Pasten et al., 2014), however, more recently accumulation models for soft soils have also been presented (Ni et al., 2014;Chen et al., 2020).The majority of the accumulation models for soft soils still lack fundamental features of natural soft soils, such as rate-dependency, anisotropy and degradation of bonding (Leroueil and Vaughan, 1990;Wheeler et al., 2003;Mitchell and Soga, 2005;Karstunen et al., 2005).
In order to overcome the limitations of the models found in the literature, a new constitutive model, Creep-Sclay1Sc, is developed to capture the effects from cyclic loading of soft soils.This is done by extending the original Creep-Sclay1S model (Sivasithamparam et al., 2015;Gras et al., 2018) by means of an explicit cumulative viscoplastic multiplier.The framework of the model is similar to the accumulation models found in literature, i.e. the combination of an existing constitutive model to simulate the first load cycle (Creep-Sclay1S) with an accumulation function (viscoplastic multiplier).With the proposed model, Creep-Sclay1Sc, the aim is to capture the change in creep rates due to long-term cyclic loading with a reasonable accuracy, and not the development of a fully-fledged model that captures all intricate details under arbitrary loading histories as function of time.
The Creep-Sclay1S model (Sivasithamparam et al., 2015;Gras et al., 2018) was chosen as a starting point of this research, as it is a viscoplastic rate-dependent model that accurately describes fundamental features of natural soft soils, such as anisotropy, structure and ratedependency (creep).Simulations of statically loaded test embankments with Creep-SClay1S model have shown that the consolidation and creep behaviour of natural soft soils is captured by the model with high fidelity (Sivasithamparam et al., 2015;Amavasai et al., 2018;Tornborg et al., 2021).
The first part of this paper presents the description of the Creep-SCLAY1Sc model and discusses the determination of the input parameters.The latter is done for the standard Creep-Sclay1S, as well as its extension for cyclic loading.The second part of the paper shows two applications of the model, within the framework of finite element analysis.The first application discusses the performance of the model by comparing the results with experimental undrained cyclic triaxial tests performed on Onsøy clay (Wichtmann et al., 2013).The second application illustrates the robustness of the model by simulating an embankment on soft soil, subjected to cyclic periodic loading, and discusses the effects that loading frequency and amplitude have on the (long-term) embankment response.

Base model: Creep-Sclay1S
The cyclic accumulation model extends the viscoplastic Creep-SClay1S model as presented in Sivasithamparam et al. (2015) and Gras et al. (2018).This model was chosen as a starting point, because it accounts for fundamental features of natural soft soil behaviour, such as anisotropy, bonding and rate-dependence.The Creep-Sclay1S model is illustrated in Fig. 1 in the p ′ -q plane.The stress state of the soil is defined by three reference surfaces: Normal Consolidation Surface (NCS), Current Stress Surface (CSS) and Intrinsic Compression Surface (ICS).The NCS is analogous to the bounding surface and delimits small and large creep strain rates, being defined by the (isotropic) pre-consolidation pressure p ′ m .CSS represents the current mean effective stress state, and is defined by the mean hydrostatic effective stress p ′ eq .ICS is introduced to represent the effect of bonding, following the framework proposed in Gens and Nova (1993), and represents the size that NCS surface would have if there were no bonding but the void ratio would be the same.Thus, ICS is an imaginary surface of which the size is defined by the intrinsic isotropic preconsolidation pressure p ′ mi .In Creep-SClay1S an associated flow rule is adopted.Creep-Sclay1S is only applicable to normally and slightly overconsolidated clays.Hence, the newly proposed accumulation model inherits this model limitation.
One of the main characteristics of the Creep-Sclay1S model is that there is no purely elastic region, which implies that creep deformations occur at all stress states.The model predicts a unique critical state, regardless of the strain rate and stress path, as the NCS rotates as a function of viscoplastic strain.Benchmark calculations with Creep-SClay1S model of laboratory and field tests performed, show that the creep behaviour of natural soft soils is correctly captured by the model (Sivasithamparam et al., 2015;Amavasai et al., 2018;Tornborg et al., 2021).

Formulation of cyclic accumulation
Similarly to the classical elasto-plasticity theory, the total strain tensor ε ij can be decomposed in the elastic ε e ij and viscoplastic ε vp ij strain components: The viscoplastic strain ε vp ij corresponds to the irrecoverable deformations within the material, and consists of two components: where ε c ij corresponds to the (static) viscoplastic creep strain, and ε cyc ij represents the additional viscoplastic strain that develops during cyclic loading.The decomposition of the strain tensor into elastic and viscoplastic components has already been successfully adopted to address the effects of cyclic loading [e.g.] (Suiker and de Borst, 2003;Niemunis et al., 2005;Karg et al., 2010).
The cyclic cumulative deformation is accounted for by extending the concept of a viscoplastic multiplier, as presented in Grimstad et al. (2010), with an additional term.The viscoplastic multiplier is always present as there is no purely elastic region.Akin to other viscoplastic models, the viscoplastic multiplier is not deduced from a consistency condition, but it is based on the overstress viscoplasticity theory (Perzyna, 1963;Perzyna, 1966).
The viscoplastic creep strain rate in the Creep-Sclay1S model is defined as: (3) where p ′ eq is the equivalent mean hydrostatic stress, Λ is the viscoplastic creep multiplier, and σ ′ ij is the effective stress tensor.The cyclic cumulative strain is incorporated by adding an additional viscoplastic multiplier Ω. Eq. 2 is hence rewritten as: The formulation implies that the viscoplastic strain resulting from the cyclic loading follows the same direction as the viscoplastic creep strain, as the associated flow rule is maintained.Experimental evidence from quasi-static monotonic trixial testing (Graham et al., 1983;Korhonen and Lojander, 1987) suggests that the assumption of an associated flow rule is a reasonable approximation for natural clays when combined with an inclined yield surface.The model is suitable to simulate a large number of cycles under quasi-static conditions; i.e. the cyclic loading period is sufficiently long to disregard the inertial effects in the soil caused by the cyclic loading.

Cyclic viscoplastic multiplier
The derivation of the cyclic viscoplastic multiplier has been performed by analysing cyclic undrained triaxial tests performed on natural soft clays in the laboratory (Zhou and Gong, 2001;Li et al., 2011;Wichtmann et al., 2013;Wang et al., 2013;Ni et al., 2014).
Fig. 2 illustrates a range of the undrained triaxial cyclic tests presented in the literature (Li et al., 2011;Wichtmann et al., 2013).In the figure q cyc represents the deviatoric cyclic stress amplitude, p ′ 0 the initial mean effective stress, S u the monotonic undrained shear strength and T the loading period.The analysis of the data from cyclic triaxial tests (Zhou and Gong, 2001;Li et al., 2011;Wichtmann et al., 2013;Wang et al., 2013;Ni et al., 2014)

yield the following observations:
• There is a semi-logarithmic relation between the development of axial strain and the number of cycles (or time), until the onset of failure; • The number of cycles required to reach failure reduces with increasing cyclic shear stress amplitude; • The number of cycles to failure increases with increasing loading frequency.
Based on these observations, the cumulative viscoplastic multiplier Ω was derived by adapting the Bailey-Norton power law (Kobelev, 2014), commonly used to characterise viscoplasticity: where c, m and σ m represent material parameters, and σ the stress (Bråthe and Josefson, 1979).
The derived cumulative viscoplastic multiplier Ω yields: where ζ, Ξ, ι and t ref|0.1% are model parameters for cyclic accumulation, T is the loading period, T 0 is the reference period, q cyc is the amplitude of the cyclic deviatoric stress and M(θ) the stress ratio at critical state, which is a function of Lode angle θ. p ′ and q represent the current mean effective stress and deviatoric stress, respectively.Eq. 6 establishes the phenomenological relation between the cyclic strain rate and the experimental observations, by means of three normalised terms.Term i defines the gradient of the cyclic strain; term ii defines the frequency dependency; and term iii defines the relative importance of the cyclic loading in relation to the static loading.
The model parameter ζ can be interpreted as the gradient of the cyclic strain measured at a reference time t ref|0.1% , at which the deviatoric strain reaches a value of 0.1%.The parameter Ξ expresses the relation between the cyclic tests performed at a loading period T and the cyclic tests performed at a reference period T 0 .From the experimental data it was clear that a loading period dependency exists, hence the need for a parameter accounting for it in the viscoplastic multiplier.The denominator pM(θ) − q in Eq. 6, represents the (deviatoric) distance of the current stress state to the NCS surface.As the distance between the CSS and NCS decreases, the accumulation rate will increase.The model introduces the notion of a reference time, which is similar to the definition of a reference creep time, as introduced in basic creep models (Leoni  Wichtmann et al., 2013Wichtmann et al., ). et al., 2008)).This poses no limitations into the applicability of the model in current practice, as shown in benchmark calculations with Creep-SClay1S model (Amavasai et al., 2018;Tornborg et al., 2021), even for a complex stress history.Combining the (apparent) preconsolidation pressure and the reference time for the laboratory test allows to scale the time-dependent response at boundary value level.

Parameter derivation for natural Onsøy clay
The Creep-SClay1Sc model contains 19 model parameters, of which 15 are shared with Creep-Sclay1S as presented in Sivasithamparam et al. (2015) and Gras et al. (2018).The parameters for Creep-SClay1Sc can be grouped into five different categories (see Tables 1 and 2 for full details on the model parameters): • Isotropic parameters related to the isotropic compressibility and the ultimate strength (λ * i , κ * , ν, M c , M e , OCR, e 0 ), similar to the Modified Cam-Clay model (Roscoe and Burland, 1968).The subscript i refers to an intrinsic value.Different values for the stress ratio at critical state in compression and extension, M c , M e , respectively, enable accounting for Lode angle dependency, see Sivasithamparam et al. (2015); • Anisotropic parameters (ω, ω d , α 0 ), similar to the S-Clay1 model (Wheeler et al., 2003) describe the evolution of anisotropy and initial anisotropy; • Destructuration parameters (a, b, ξ 0 ), similar to the S-Clay1S model (Karstunen et al., 2005) describe the degradation of bonding and the initial amount of bonding; • Viscosity parameters (μ * i , τ), similar to the Creep-SClay1S model (Sivasithamparam et al., 2015;Gras et al., 2018) describe the intrinsic creep rate and the reference time used in the test for determination of OCR; • Cyclic accumulation parameters (ζ,Ξ,ι,t ref|0.1% ), as presented in Eq. 6.
All necessary parameters for the description of Creep-Sclay1Sc, with the exception of the cyclic parameters, can be derived from standard geotechnical laboratory tests, such as the combination of undrained monotonic triaxial tests and oedometer tests.The additional parameters related to the cyclic accumulation require the execution of additional undrained cyclic triaxial tests.Although, the number of model parameters is significant, the Creep-SClay1S has been successfully adopted in industry applications (POVM, 2020;Tornborg et al., 2021).In the current analysis, the results on natural Onsøy clay were used to evaluate the performance of the cyclic model, due to the considerable amount of experimental data available from high quality block samples of natural clay with extensive static and cyclic test programmes (Andersen, 2009;Berre, 2014;Grimstad et al., 2014;Wichtmann et al., 2013).
The isotropic, anisotropic and destructuration parameters were determined from monotonic triaxial and oedometer tests performed on the natural Onsøy clay (Jostad and Berre, 2010).The viscosity parameters μ * i and τ, were determined from triaxial tests on natural Onsøy clay, as documented in Jostad and Yannie (2015).Table 1 presents a summary of the adopted parameters.The derivation of the cyclic parameters will be addressed in Section 4.1.

Monotonic element tests on natural Onsøy clay
Fig. 3 shows the comparison between the experimental data and the numerical results (at a single Gauss point) for monotonic triaxial tests (Wichtmann et al., 2013) and constant rate of strain compression tests CRSC (Jostad and Berre, 2010).The numerical results are in reasonable accordance with the experimental results.Fig. 3a shows that for the triaxial tests in compression, the peak shear stress and slope of the critical state line are predicted well.Similarly, the response in triaxial extension are reasonably well captured.In both tests, the samples have been initially anisotropically consolidated to the estimated in situ stress level.The experimental results indicate that the stress path is not completely undrained, given there is a reduction of the mean effective stress p ′ during the initial loading stage.This might be related to either sample disturbance, the compressibility of the soil skeleton, or most likely, indicates that the sample has not been fully saturated.In the triaxial extension path, the anisotropy changes considerably, affecting the predicted effective stress path, as discussed for the SClay1 model in Wheeler et al. (2003).Fig. 3b illustrates that for a one-dimensional compression tests both elastic and viscoplastic components of the axial strain, as well as the preconsolidation pressure, are accurately predicted.In addition, the rate of destructuration is captured reasonably well.

Determination of cyclic accumulation parameters
The cyclic accumulation parameters (ζ, Ξ, ι and t ref|0.1% in Eq. 6), were derived from undrained cyclic triaxial tests performed on natural Onsøy clay (Wichtmann et al., 2013).The reference time t ref|0.1% was found to be related to the ratio of the deviatoric cyclic stress amplitude q cyc and the initial mean effective stress p ′ 0 (Fig. 4).Assuming a power relation between t ref|0.1% and q cyc /p ′ 0 , Eq. 6 is reformulated as: where Γ α and Γ β result from a non-linear least squares fit of the experimental data, as shown in Fig. 4. The occurrence of negative values for the reference time are prevented by using a power equation, which is a desirable feature from a computational point of view.The mean square error of the fit to the experimental data is 5.53 × 10 − 4 .The model parameter Ξ captures the rate effect of the loading frequency (period) on the cyclic accumulation.Fig. 5a shows the strain rate for three cyclic triaxial tests performed under similar initial effective stress and loading conditions, at different loading periods.The results suggest that a relation between the loading period and the deviatoric strain rate can be established.Fig. 5b shows the average normalised strain rate for the three loading periods above, together with the numerical results after parameter optimisation.The average normalised strain is defined as the average ratio between the strain at the loading period T and at the reference period T 0 , until the onset of failure.
Although the experimental and numerical results agree well, these results are based on the analysis of only three cyclic tests.Therefore, more tests should be conducted in order to further investigate the effect of the loading period.
The model parameters ζ and ι were obtained by means of monoobjective optimisation of the cyclic triaxial tests (Gras et al., 2017).The values adopted for the analyses in this paper only considered the triaxial tests that recorded more than 30 cycles to failure.ι turns out to be one order of magnitude smaller than parameter ζ, with a similar amount of scatter.This scatter is most probably related to the fact that the Onsøy clay is a natural clay.Substantially more test data would be required to identify clear trends in the data that are not affected by the typical scatter induced by the natural variation, sampling and testing.

Cyclic triaxial tests on natural Onsøy clay
The performance of the Creep-SClay1Sc in modelling of cyclic loading was assessed by comparing the numerical results with experimental undrained cyclic triaxial tests (Wichtmann et al., 2013), at a single Gauss point level.Fig. 7 presents the comparison between the experimental data and the numerical results for four different cyclic triaxial tests.The axial strain corresponds to the permanent axial strain developed during the cyclic loading.The model parameters adopted for the analyses have been presented in Tables 1 and 2 for the Creep-Sclay1Sc model.Table 3 summarizes the description of the cyclic triaxial tests in the laboratory that are used for comparison.
In general, the numerical and experimental results compare favourably.Fig. 7c illustrate that the new model is able to correctly predict the accumulation of axial strain in a cyclic compression test until the onset of failure.The moment the failure occurs, the experimental results show that the strain develops at a higher rate than predicted by the numerical simulation.accumulation of axial strain, as well as the onset of failure are again correctly captured by the Creep-SClay1Sc model.Therefore, the rate effect resulting from the loading period is properly incorporated in the model.
For test 23 (Fig. 7c), which was loaded at isotropic average stresses, no failure was recorded or predicted.The isotropic consolidation paths in the latter test gradually destroys all bonds in the natural samples, and when combined with evolving anisotropy, it is likely that the latter largely cancels out the generation of additional pore water pressures resulting from the gradual collapse of the fabric.
Similarly to the monotonic test results, the simulation of the triaxial extension test was not as accurate as the predictions for compression tests (Fig. 7d).Both the strain rate as well as the number of cycles required to the onset of failure are overestimated by the Creep-SClay1Sc model.Since data from only one triaxial test in extension was available, it is not possible to assess whether the difference between the numerical and experimental results are representative.There could have been issues, for example, due to inhomogeneities within the sample, the necking of the sample, membrane effects, or other unknown experimental inaccuracies.

Numerical implementation
The Creep-Sclay1Sc model is an explicit model (Niemunis et al., 2005), in the sense that the accumulative plastic strain is not computed for each load cycle but is instead modelled by means of the viscoplastic multiplier.The use of explicit models at boundary value problems, require a two-step approach (Niemunis et al., 2005;Pasten et al., 2014), in which initially one reference load cycle is simulated to compute the cyclic stress amplitude, followed by the application of the cumulative model (see Fig. 8).The procedure for the execution of boundary value analysis is as follows: 1. Calculation of initial stress field (K 0 procedure), caused by the self weight with Creep-Sclay1S followed by the initialisation of the state variables in each integration point; 2. Calculation of any monotonic loading history, e.g.embankment construction followed by consolidation, with Creep-Sclay1S; 3. Calculation of one reference load cycle with Creep-Sclay1S and computation of the cyclic deviatoric stress amplitude, q cyc , at each integration point; 4. Calculation of stress and strain with the accumulation model, Creep-Sclay1Sc, at each integration point, assuming that cyclic deviatoric stress amplitude remains the same during the calculation.If required this step can be interrupted, and a new reference load cycle can be applied to simulate a subsequent change in the cyclic loading conditions.
By following the procedure above, the equilibrium is guaranteed during the entire calculation, according to the adopted boundary conditions for the analysis.

Cyclic triaxial test
The triaxial tests presented in Section 4.2 have again been modelled, but this time at boundary value level, using the finite element package Tochnog Professional (2017).A 2D axisymmetric modelling approach was used to model the undrained cyclic triaxial test on a sample of 54 mm diameter and 108 mm height.Fig. 9 shows the dimensions and the mesh for the triaxial test simulations.The domain is meshed with a total    of 217 quadrilateral elements using second order shape functions, with 8 nodes in each element.The displacements were constrained along the horizontal direction along the centre axis, and vertically along the bottom edge.The top and outer boundary have a distributed load that is in equilibrium with the anisotropic stress state after consolidation.All flow boundaries are impermeable, and the initial pore pressures are set at 600 kPa, following the experimental description presented in (Wichtmann et al., 2013).
The model parameters for the material are the same as presented in Tables 1 and 2. The analyses have been done using a coupled formulation between the solid skeleton and groundwater.The adopted groundwater parameters are the following: water density of 1000 kg/ m 3 ; bulk stiffness for the water of 10 6 kPa; and hydraulic conductivity of 10 − 9 m/s.The remaining cyclic loading conditions correspond to test 15, which is already described in Fig. 7a and Table 3.
The result of the boundary value analysis is shown in Fig. 10 with the experimental data and the previously shown single point results as comparison.The initialisation of the cyclic component of the load is visible, as well as the correctness of the modelling strategy.The small differences between the single point and boundary value simulations are expected, given absence of localisation at small strain magnitudes that is further promoted by the regular mesh and the axisymmetric boundary conditions used.At larger strains there is a discrepancy between the numerical and experimental results.These can be explained by the effects of strain localisation, and the formulation of the finite element method, which is not well suited to address problems involving large deformations.Furthermore, the experimental data used in the comparisons also is less reliable after strain localisation occurs as it leads to a non-uniform stress state that precludes boundary level stress and porewater pressure measurements.
After the successful verification of the model implementation at boundary value level using a fully coupled formulation the model is ready to be used in a more demanding boundary value problem.

Embankment on soft soil
In order to illustrate the benefits of using the proposed cyclic accumulation model for modelling the effect of a large number of load cycles at boundary value level, a generic embankment on Onsøy clay has been simulated.Although the focus of the analysis is to illustrate the performance and applicability of the new model, the geometry and loading conditions of the problem resemble a simplification of gravity foundations for onshore wind turbines, machine foundations, and road and railway embankments under traffic loading conditions.
The embankment has a height of 2 m, a crest width of 12 m and slope inclination of 1:2, and it is assumed to rest on top of a dry crust of 4 m thickness.The natural subsoil consist of a natural soft clay layer with a thickness of 36 m.The embankment and dry crust are assumed to behave as linear elastic materials (material properties in Table 4), whilst the clay parameters adopted correspond to the Onsøy clay, and follow the newly proposed model (material properties in Tables 1 and 2).
Half of the embankment problem was simulated under plane strain conditions using an Updated Lagrangian formulation, by utilising a symmetry condition.The domain was discretised using 2 nd order triangular elements (in total 11 095 elements).The groundwater level was assumed to be placed at the interface between the dry crust and the soft clay layer below.The horizontal displacements are constrained at the two side boundaries, whereas the vertical and horizontal displacements are constrained at the bottom boundary.The hydraulic boundaries are closed at the sides and bottom with a free water surface just below the dry crust.A complete overview of the geometry and mesh of the simulated embankment is shown in Fig. 11.
The system of equations was solved following a fully coupled   consolidation formulation, in which the groundwater flow is accounted for by means of the storage equation.The analysis was performed using the finite element package Tochnog Professional (2017).
The cyclic loading was simulated by applying a distributed harmonic load (frequency 1 Hz) slightly out of centre on top of the embankment, with a 2 m width and 500 kN/m amplitude.Following the procedure explained in Section 5.1, the load was applied first as one cycle, to capture the amplitude of the cyclic deviatoric stress throughout the domain, which was followed by the application of the cyclic accumulation (as described in Fig. 8).
Fig. 12 shows the evolution of the vertical displacements as function of the number of cycles for the point located at the interface between the embankment and the dry crust, along the axis of symmetry.The results of the basic creep model, Creep-SClay1S, without cyclic accumulation are also plotted for comparison.The results suggest that the cyclic loading has a significant effect on the vertical displacement.For up to 10,000 loading cycles the responses predicted by the basic creep model and the new cyclic model are similar.However as the number of cycles increases (>10,000), the effect of the cumulative straining due to cyclic degradation becomes dominant.suggest that up to 10,000 load cycles there is little deformation of the embankment and subsoil (Figs.13a and 13b).As the number of loading cycles increases, the displacements increase, as well as the shear strain levels.The largest displacements occur underneath the embankment, while the largest shear strain is registered in the soft clay layer.As the  number of cycles increases, there is a failure of the embankment.This is clear from the total displacement and shear strain fields at 50,000 cycles (Figs.13e and 13f) which exhibits a classic shear failure of a foundation.The effect of the loading amplitude and the loading frequency were assessed in more detail by performing two additional computations, where either the loading amplitude, or the loading frequency have been doubled: • Case 1: Loading frequency 1 Hz, loading amplitude 500 kN/m; • Case 2: Loading frequency 1 Hz, loading amplitude 1000 kN/m; • Case 3: Loading frequency 2 Hz, loading amplitude 500 kN/m.Fig. 14 presents the development of the vertical displacement for the point located at the interface between the embankment and the dry crust, along the axis of symmetry, for the three different cases.The increase in the load amplitude (case 2) from 500 kN/m to 1000 kN/m is predicted to result in larger displacements, whereas the increase in the loading frequency (case 3) leads to a smaller displacement, at any given number of cycles.This is consistent with the model formulation and with the experimental findings reported in Li et al. (2011) andWichtmann et al. (2013).
The simulations with the cyclic accumulation model suggest that doubling of frequency has a bigger impact on the predicted results than the doubling of the load.This is because, for this particular application, when the load is doubled, the increase in the cyclic deviatoric stress on the soft soil layer is relatively small, when compared with the initial at rest mean and deviatoric stresses (due to the existence of the embankment and dry crust).These results may however differ, for other types of applications.
It should be noted that the predicted displacements in these analyses are relatively high.This is related to the low strength and stiffness of the Onsøy clay tested (both site conditions and block samples from shallow depth), which is assumed to apply for the complete clay layer, without correcting for the in situ OCR profile.In Fig. 15 the vertical displacement for several points located at the interface between the embankment and dry crust layer (Fig. 15a), and at various depths along the vertical axis of symmetry (Fig. 15b) are plotted.A gradual reduction of the vertical displacement at the interface between the embankment and the dry crust is found in all the analyses, from the axis of symmetry towards the edge of the embankment.Again, Case 1 is found to lead to larger displacements than Case 2. Up to 50,000 cycles the displacements of Case 3 are nearly negligible.
Regarding the variation with depth (Fig. 15b), all the analyses show that the points located in the dry crust (levels between 36 and 40 m) have a smaller attenuation with depth, whereas the attenuation significantly increases for the points located in the clay layer (level below 36 m).This is related to the fact that the dry crust is modelled as a linear elastic material, in which no (visco) plastic deformation occurs.
Fig. 15 shows again that the doubling of frequency has a bigger impact than the doubling of the load amplitude.For 50,000 cycles the maximum vertical displacement is 0.08 m for case 3, whilst this value reaches 0.75 and 0.80 m, for cases 1 and 2, respectively.The difference in the predicted vertical displacement for Cases 1 and 2 tends to reduce with the increase in depth (Fig. 15b).This is due to the spreading of the load: as the depth increases, the relative deviatoric stress increment from the cyclic load also decreases.

Conclusions
The paper proposes a new constitutive model developed for cyclic accumulation.The proposed model is able to implicitly capture the effects from cyclic loading of soft natural soils, by extending the original Creep-Sclay1S model.The strain accumulation due to cyclic loading is explicitly accounted for by means of a viscoplastic multiplier, in a similar manner to the viscoplastic effects associated to rate-dependence in the original model.The model extension adds four additional model parameters, which can be determined from undrained cyclic triaxial test data.The model accounts for the effect of the cyclic loading amplitude, as well as the loading period i.e. loading frequency.This approach leads to a computational efficient and numerical stable solution for the simulation of a large number of load cycles without encountering significant numerical issues.
The comparisons of the model simulations and experimental data from undrained cyclic triaxial tests on block samples of natural Onsøy clay demonstrate that the model is able to predict accurately the rate of strain accumulation, as well as the onset of failure, for isotropically and anisotropically consolidated triaxial tests sheared in compression.For loading paths in extension the model underpredicts both the strain rate and the number of cycles to failure, i.e. the results are non-conservative.However, only one triaxial extension test was available.A power relation between the loading period and the reference loading period T ref was derived using only three cyclic triaxial tests.The parameterisation of the cyclic model exhibits some scattering.In order to improve the model reliability, and fully validate and refine the parameter determination, more cyclic testing on natural soft soil is required.These tests should be performed at different initial stresses, and along different stress paths, for a range of loading rates, both in compression and extension.It should be highlighted that a natural scatter will always be present when dealing with intact samples of natural clay that stems from the natural variability in the soil, and uncertainties induced by the sampling, transport, storage and testing procedures.
The proposed model has been shown to work on both single element level and boundary value level.The effects of the cyclic loading have been illustrated for a generic embankment built on top of a deposit of Onsøy clay.The analyses have highlighted the effect of cyclic loading on the evolving rate-dependent deformations in the soil below the embankment.It was shown that, for the case studied, the cyclic component is more significant than the static consolidation and creep.Furthermore, there is a relation between the loading amplitude and/or the loading frequency and the accumulated deformation in the subsoil.
The newly proposed model offers novel insight into the cyclic response of soft natural soils.In addition to its capability to correctly model soft soils under monotonic loading, the model offers the possibility to simulate large numbers of load cycles in an efficient and accurate manner.This enables the simulation of the long-term response of cyclically loaded foundations on soft soils, for which the serviceability limit state considerations (over a long period of time) are governing the design.Examples of the latter include machine foundations, road and railway embankments and foundations for wind turbines.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 2 .
Fig. 2. Experimental results from undrained cyclic triaxial tests.Development of axial cyclic strain: (a) at different cyclic shear stress amplitudes (adapted from Li et al., 2011), (b) at different cyclic shear stress amplitudes and (c) at different loading periods (adapted from Wichtmann et al., 2013).

Fig. 3 .
Fig. 3. Comparison of experimental and numerical results: (a) anisotropically consolidated undrained monotonic triaxial tests and (b) constant rate of strain consolidation test.

Fig. 6
Fig. 6 summarises the values of the model parameters ζ and ι for all the analysed tests.It was found that the values of ζ, with the exception of three tests, namely 16, 17 & 19, fall within the same range.The parameter

Fig. 8 .
Fig. 8. Schematic representation of the calculation procedure for the accumulation model.

Fig. 9 .
Fig. 9. Finite element mesh of the triaxial test as a boundary value problem (Dimensions are in mm).

Fig. 12 .
Fig. 12. Vertical displacement of the point located along the vertical axis of symmetry at the interface between the embankment and the soil.
Fig. 13 presents the spatial distribution of the total displacements (figures (a), (c) and (e)) and shear strains (Figures (b), (d) and (f)) after a different number of loading cycles in the full domain.The simulations

Fig. 14 .
Fig. 14.Vertical displacement of the point located along the vertical axis of symmetry at the interface between the embankment and the soil, for three different loading conditions.

Fig. 15 .
Fig. 15.Vertical displacement at different number of cycles (10,000, 10,000 and 50,000): (a) at the interface between the embankment and subsoil and (b) along the vertical axis of symmetry.

Table 1
Model parameters for the natural Onsøy clay.

Table 3
Description of cyclic triaxial tests.

Table 4
Elastic parameters for the embankment and the dry crust.