Investigation of three sophisticated constitutive soil models: From numerical formulations to element tests and the analysis of vibratory pile driving tests

The performance of three advanced constitutive models has been evaluated based on element tests and on a comparative study on the simulation of vibratory pile driving tests in saturated sand. The inspected constitutive models are the Sanisand model and Hypoplasticity with Intergranular Strain (Hypo+IGS) as well as with Intergranular Strain Anisotropy (Hypo+ISA) extension. The performance of the constitutive models is first evaluated by the simulation of element tests used for the parameter calibration of the sand used in the model tests. The constitutive models are then applied for the simulation of a vibratory pile driving test. The pile penetration, the driving force, the pore water pressure development and the incremental displacement in the vicinity of the pile tip are compared to the measurements in the model tests. The strengths and weaknesses of the different constitutive models are assessed. Generally, the model predictions showed good agreement with the experimental results. Despite different constitutive formulations (hypoplastic vs. elasto-plastic), all three models were able to reproduce the main mechanisms of the driving process properly. It may be concluded that all three models allow a proper prediction of vibratory pile driving as long as a proper calibration of the material parameters is secured.


Introduction
Numerical simulations of geotechnical structures often require the use of sophisticated constitutive soil models, especially if the investigated structure is subjected to cyclic loading. These constitutive models are expected to reproduce the mechanical behaviour of the soil at different stress/strain amplitudes, particularly under undrained cyclic loading, whereby cyclic mobility effects and the accumulation of pore water pressure should be adequately described. In order to accurately address these effects the complexity of advanced constitutive models grows exponentially. Beside the well-known pyknotropy (void ratio dependence of the soil behaviour), barotropy (pressure dependence), viscosity (time-dependence) and anisotropy of some types of soils, the main challenge has turned out to be the reproduction of the specific soil behaviour under cyclic loading. For this purpose different mechanisms have been introduced into the different families of constitutive models, such as the history variables/surfaces and their evolution laws incorporated into elastoplastic frameworks (Dafalias and Manzari, 2004;Benz, 2007;Tafili, 2020) or the concept of the Intergranular Strain into the Hypoplasticity (Niemunis and Herle, 1997;Fuentes and Triantafyllidis, 2015), amongst others. Today many different constitutive models coexist and modifications/extensions as well as novel formulations are proposed regularly. However, due to the complexity of the models, performance checks are often only conducted based on element tests, see e.g. (von Wolffersdorff, 1996;Niemunis and Herle, 1997;Herle, 1997;Dafalias and Manzari, 2004;Mašín, 2013;Fuentes et al., 2017;Bode et al., 2020;Tafili and Triantafyllidis, 2020a,b). Likewise, the direct comparison of those soil models is typically also performed based on element tests, e.g. (Wichtmann et al., 2019;Tafili et al., 2020). A comparison of the prediction quality in boundary value problems is often missing but attracting growing interest, see e.g. (Jostad et al., 2020;Roy et al., 2020). The missing assessment is not surprising, as the simulation of boundary value problems demands a very thorough implementation of the soil models to assure sufficient numerical stability and robustness. To ensure the required robustness, often (undocumented/non-communicated) slight modifications of the constitutive model and detailed expertise are required. In addition, these advanced constitutive models are often associated with time-consuming and complex parameter calibration procedures requiring a large number of tailored laboratory tests. One of this paper's objectives is to demonstrate that a parameter calibration for the models at hand can be based on a few well chosen laboratory tests to reach a satisfactory prediction of advanced boundary value problems.
Since vibratory pile driving in saturated soil involves cyclic loading, possible liquefaction and large volumetric changes in the soil, the application of different constitutive models to such a boundary value problem is able to provide a variety of comparisons and avoids the risk to proof the suitability for one specific topic while being unable to capture other relevant aspects. So far, the numerical simulation of vibratory pile driving in water-saturated soil has been subject of only a few studies (Galavi et al., 2017;Chrisopoulos and Vogelsang, 2019;Staubach and Machaček, 2019;Staubach et al., 2020b;Giridharan et al., 2020). An adequate numerical modelling approach for pile driving is not only relevant for the prediction of the pile penetration and soil resistance during the installation process but also for the study of installation-induced changes in the state of the soil surrounding the pile. Different studies showed that the installation process can have a distinctive influence on the pile response when subjected to vertical or horizontal loading following its installation (Heins and Grabe, 2017;Fan et al., 2019;Staubach et al., 2020a) and that its incorporation in the numerical model is mandatory in order to predict the response of the pile correctly (Fan et al., 2021;. For the simulation of the vibratory pile driving model test considered in this work, three sophisticated constitutive models are compared: Hypoplasticity with Intergranular Strain (Hypo+IGS) extension and with Intergranular Strain Anisotropy (Hypo+ISA) extension as well as the Sanisand model. Due to the large deformations during the pile penetration, an updated Lagrangian formulation is employed.
The paper is structured as follows: First, the vibratory pile driving model tests performed by Vogelsang (2017) are introduced in Section 2. The corresponding numerical model is described in Section 3 and complemented by details of important aspects of the numerical model in Section 4. The three constitutive models used in this study are introduced in Section 5. This also includes notes on modifications of the governing equations and implementation strategies which deviate from their original publication. The parameter calibration is subject of Section 6, while Section 7 evaluates the performance of the soil models based on element tests. Section 8 contains the comparison of the simulation results with the measurements from the model tests. The paper closes with a summary and a conclusion in Section 9.

Model tests
A detailed description of the vibratory pile driving experiments is given in Vogelsang (2017). The most important aspects will be recapped here to disclose the general concepts. The experiment discussed in this work was carried out in half of an axisymmetric test bench with a diameter of 0.94 m and an acrylic glass plate at the front (a schematic illustration as well as a picture of the device is given in Fig. 1). The actual height of the embedded soil varies only little and is idealised as constant h Sand = 0.81 m in the recalculations of the tests. The pile has a diameter of d Pile = 33 mm and a 60 • pointed tip.
Throughout the experiments the so-called "Karlsruhe Sand" has been used with its index parameters given in Table 1. It is important to note that during the last decades, different sands have been referred to as "Karlsruhe Sand". An extensive experimental study and description of the "current Karlsruhe Sand" is given in Vogelsang (2017). These experiments form the basis of the calibration of the material parameters for the back calculations. The Kozeny/Carman-relation linking the hydraulic conductivity k to the porosity n of the material reads (Kozeny, 1927;Carman, 1939): Therein, γ w = 10 kN/m 3 is the specific weight and η w the dynamic viscosity of the pore-water. The effective grain size of "Karlsruhe Sand" is d e = 0.5 mm and C = 308 (Vogelsang, 2017). As indicated in Fig. 1(b), the pile was installed after the soil was placed. The dry sand was pluviated into the model container filled with deaerated water. Hammer impacts on the container base were used to compact the sand to its target density, the "initial density" of the pile installation tests. The height of the groundwater level corresponds to the top of the soil. Vogelsang (2017) states that this method results in nearly saturated, homogeneous samples. The initial relative density of the present model test was D r0 = 71%, which corresponds to an initial void ratio of e 0 = 0.637. Then, the pile was monotonically pushed into the Fig. 1. Schematic illustration of the test device (a) and picture of the device (b) with the already embedded soil. A view from above is given in figure (c). The pore pressure transducers (PPT) are marked according to their location (not scaled, based on Vogelsang (2017)). a The critical friction angle was determined as the inclination of a loosely pluviated cone of sand, i.e. as the angle of repose. soil to its "initial depth" of 0.16 m. Starting from this configuration, the driving force with a frequency of f = 25 Hz was applied using a vibrator mounted guide-free on top of the pile. The static moment of the vibrator was M stat = 5.33⋅10 −3 kg⋅m and the combined vibrating mass of the pile, the load cell and the vibrator is m vib = 7.881 kg.
The measurement concept includes the recording of the pile head displacement and the pile head force measured by means of a load cell installed between the pile head and the vibrator. Changes in pore pressure were measured at two pore pressure transducers PPT A and PPT B installed at the front glass plate as displayed in Fig. 1. Referred to the driving process in the experiment, PPT A is passed by the pile tip after a few seconds of vibration whereas PPT B is reached after approximately 6 s. In addition, displacement fields in the soil were obtained from digital image correlation (DIC) analyses.

Numerical model
The finite element simulations were performed using the finite element code numgeo 1 . The finite element mesh as well as information on the boundary conditions and dimensions are given in Fig. 2.
Exploiting the symmetry of the model test, an axisymmetric FE model was used. The forces and masses were scaled to fit to the halfaxisymmetric experimental set-up. The soil was discretised using approximately 2500 quadratically interpolated elements, resulting in a mean nodal distance of about 15 mm near the symmetry axis. The soil was discretised using so-called u-p elements, with the soil displacements (u) and the pore water pressure (p) as primary unknowns. The so-called zipper-method was used to avoid mesh distortion when the pile penetrates into the soil. This method is well established for such axisymmetric penetration problems (van den Berg, 1994;Henke and Grabe, 2008). By using this approach, the (left) boundary of the soil in the symmetry axis below the pile tip is not constrained in horizontal direction by Dirichlet boundary conditions, but by a contact constraint with a thin vertical extension of the pile directly in the symmetry axis (see the red line in Fig. 2). This allows the finite elements below the tip to be pushed to the side when the pile penetrates into the soil. In addition, the displacements of the soil at the right boundary are constraint in horizontal direction and at the bottom boundary in vertical direction. As indicated in Fig. 2, the pore water pressure at the top surface of the soil is prescribed (p w = 0 kPa) and constant throughout the simulation. Between the pile and the unbalances, which are modelled as mass points, a spring element with the same stiffness as the load cell depicted in Fig. 1 is used to obtain the force between pile and unbalances.
A frictional contact between the aluminium pile and the soil is considered using a Coulomb friction model (see Section D). Note that instead of a perfectly pointed pile tip, as used in the experiments, the wire is rounded in the numerical model to prevent stress concentration, see Fig. 2.
Accounting for minor air inclusions estimated to < 0.1% of the water volume (degree of saturation S⩾99, 9%) originating from the sand pluviation, the bulk modulus of the pore water is assumed to be K w = 1.2⋅ 10 4 kPa. Due to the small scale of the model tests and thus small stresses in the vicinity of the pile, grain crushing has not to be accounted for.
In the numerical simulation, the initial configuration corresponds to the geostatic equilibrium with a lithostatic stress state (the soil placement process was not considered). The pile is modelled in its initial penetration depth of approximately 0.16 m (prior to the vibratory pile driving). In the experiment, the pile was monotonically pushed to this depth, which is not taken into account in the present simulations. In Staubach et al. (2020b), the installation of the pile to the depth 0.16 m was numerically simulated, showing no significant influence on the subsequent vibratory pile driving.
The simulation of the model test is conceptually divided into four consecutive steps: 1. First, the initial stress state due to the self-weight of the soil and the pile is applied in a so-called geostatic step, without generating any deformations. During this step, all displacements of the pile are constrained. For the calculation of the initial stresses, an earthpressure coefficient of K 0 = 1 −sin(φ) was assumed. For the hypoplastic models, the initial values of the Intergranular Strain tensor are set to 5⋅10 −5 in the vertical direction, thus corresponding to the parameter R, and to zero in case of all other components. In case of Sanisand, the back-stress tensor α is initialised as the ratio of initial  numgeo (see www.numgeo.de and (Machaček, 2020;) is a finite-element program, developed by the first two authors with the emphasis on the solution of non-linear, coupled (dynamic) geotechnical boundary value problems. deviatoric stress tensor and mean effective stress (the initial state is inside the elastic locus). An artificial compressive traction t stab = 0.1 kPa is applied on the top surface for numerical purposes (note that the magnitude remains constant during the following steps). 2. Then, the pile is released, activating the contact between the pile and the soil. 3. In the third step, the vibratory pile driving is simulated. The vibrator force over time and corresponding frequency are illustrated in Fig. 3. As the imbalances required a time to reach the desired frequency in the experiments, a linear increase in the force amplitude was assumed for the first 0.075 s. A constant time increment of Δt = 0.001 s and a numerical damping for the Hilber-Hughes-Taylor (HHT) time integration scheme of α = −0.05 are used during the dynamic step. The step time is 6.35 s, which equals the time of vibration considered in the experiment. 4. In order to take into account the resting phase following the driving, a second dynamic step with a step time of 0.65 s is added. The time integration scheme and increment size are the same as in the previous step.

Soil as a two-phase porous medium
For this work, the soil is idealised as a two-phase porous medium consisting of grains and pore water. For the simulations the finite element code numgeo has been used, which offers various element formulations based on the Theory of Porous Media (TPM) for the simulation of two-phase and three-phase porous media in transient and dynamic analyses. In a preceding work Staubach and Machaček (2019) investigated the applicability of the u-p formulation (the arising linear system is sought for the solid displacements u and the pore water pressure p w (Zienkiewicz et al., 1980)) to the simulation of vibratory pile driving in sand. It was shown by means of a semi-analytical solution as well as a comparative numerical study that the relative acceleration between the solid phase and the fluid phase can be neglected in the present study (which is the case when applying the u-p formulation). The governing equations used for the description of the two-phase porous medium consist of the balances of linear momentum of the overall mixture, the balance of linear momentum of the pore water and balances of mass of the solid and pore water, respectively. For a detailed derivation of the implemented equations, the interested reader is referred to Staubach and Machaček (2019), where the application of the u-U and u-p-U element formulations to the simulation of the vibratory pile driving tests is presented as well.

Updated Lagrangian formulation
Due to the large deformations during the pile penetration, an updated Lagrangian formulation is employed. For this purpose the Jaumann-Zaremba stress rate is used where the objective stress rate σ was introduced. W is the spin-tensor.
To integrate Eq. (2), the Hughes-Winget algorithm (Hughes and Winget, 1980) is applied. It is acknowledged that no finite strain is workconjugated to the Jaumann-Zaremba stress rate. The strain rate ε is calculated as the symmetric portion of the velocity gradient. The rateadditive plasticity formulation is adopted (Nemat-Nasser, 1982) in case of the elasto-plastic models, which is suitable for materials with small elastic range such as soils (Fish and Shek, 2000).

Linear viscosity
For very fast deformations at vanishing mean effective stress (e.g. p = |tr(σ)/3| < 2 kPa), a small viscous stress σ vis is added to the constitutive stress σ c : The viscous stress is calculated from σ vis = λ1tr(ε) + 2με. To ensure a smooth transition, λ and μ are assumed to increase linearly with vanishing mean effective stress p: where λ min = 0.1 kPa⋅s, λ max = 1.0 kPa⋅s, μ min = 0.1 kPa⋅s and μ max = 1.0 kPa⋅s are material parameters. Up to now, no experiments have been conducted to calibrate the material parameters for granular soils during the phase transition into a suspension. Therefore, the viscous stress in this work has to be regarded as an artificial stress stabilising the stressstrain behaviour of granular soils at small stress states and high deformation rates.

Constitutive models
The three constitutive models considered in the simulations and the calibration of their parameters are discussed in the following.

Hypoplastic model with Intergranular Strain
The basic hypoplastic model of von Wolffersdorff (1996) interrelates the objective effective stress rate σwith the strain rate ε: L is a fourth order tensor being linear in ε, whereas N is a second order tensor being non-linear in ε. Both stiffness tensors L and N are functions of effective stress and void ratio.
For a detailed presentation of the equations for L and N it is referred to von Wolffersdorff (1996). To improve the performance of the hypoplastic model in the range of small strain cycles Niemunis & Herle (Niemunis and Herle, 1997) introduced a new tensorial state variable, the Intergranular Strain h, which memorises the recent deformation history. With this extension the basic equation of the constitutive model reads: More details on the constitutive relations of this extension can be found in Niemunis and Herle (1997). The present implementation uses an adaptive explicit Euler scheme to integrate the stress rate within a substepping method and error control (see the modified Euler method used in Ding et al. (2007) proposed by Sloan (1987)). The error of the explicit scheme within every subincrement is calculated and the substep size is reduced if the error is too large. Likewise, the substep size is increased if the error is small.
As mentioned in the introduction some modifications to the original constitutive equations have been adopted to improve the robustness of the constitutive model. These modifications are briefly summarised in the following: 1. For tensile mean effective stress, the material behaviour is not defined. Therefore, the following correction is applied for effective stresses smaller that σ HP = 0.01 kPa: a hydrostatic stress Δσ HP ii is added to the current effective stress σ HP to guarantee (σ HP ii + Δσ HP ii )/3⩾σ HP . This modification was also applied in Chrisopoulos et al. (2016). 2. For void ratios e exceeding the highest possible void ratio e i (p) corresponding to the loosest state at a given mean effective stress p, a special treatment originally developed by A. Niemunis is applied: Therein, f ei is an additional scalar factor applied to Eq. (A.6) and Eq. (A.7) in Appendix A.

Sanisand
The constitutive model family evolving around the Sanisand version of 2004 (Dafalias and Manzari, 2004) is constantly growing in number. Many extensions have been proposed. See e.g. Lashkari (2010) for an extension with anisotropic elasticity, Taiebat and Dafalias (2008) for a yield surface with closed cap and Dafalias and Taiebat (2016) for an extension without elastic range. For an extension by a memory surface to allow for the simulation of several thousands of loading cycles the reader is referred to  and for an enhanced performance under cyclic loading in general to (Liu et al., 2020). In addition an extension by a memory surface and semifluidized states was proposed in (Barrero et al., 2020;Yang et al., 2020) and Petalas et al. (2020) taking into account the evolving fabric anisotropy of sand. Despite these extensions, the original version remains a very frequently used constitutive model due to its proven robustness in various numerical studies (see e.g. Kementzetzidis et al. (2019) and Staubach et al. (submitted for publication)). Sanisand is an elasto-plastic model of the following form: Details on the state variables α and z as well as on the equations defining the elasto-plastic stiffness E ep are provided in Appendix B. No modifications of the constitutive equations as proposed in Dafalias and Manzari (2004) have been adopted.
Different implementations exist and are partly freely available (e.g. the implementation by M. Martinelli, C. Miriano and C. Tamagnini via soilmodels.com). These implementations were tested in the course of the present work but were found to be numerically not stable for the simulation of the vibratory pile driving tests. Therefore, a new implementation has been written by the authors. Similar to the implementation of the Hypo+IGS model, an explicit substepping scheme has been used. For trial stress states outside the elastic region of the model, a return-mapping algorithm as proposed in Dafalias and Manzari (2004) was utilised. The same correction of small mean effective stresses as for the Hypo+IGS has been applied. Furthermore, an additional viscous stress rate as introduced in Section 4.3 has been adopted.

Hypoplasticity with intergranular strain anisotropy
Hereby the hypoplastic model of von Wolffersdorff (1996) is coupled with the Intergranular Strain anisotropy (ISA) (Fuentes and Triantafyllidis, 2015) as described in Poblete et al. (2016). The elastoplastic formulation of ISA enables the simulation of small strain effects by increasing the stiffness and reducing the plastic strain rate yielding the following mechanical relation between the stress rate and the strain rates: with E and ε p called the stiffness tensor and plastic strain rate under mobilised conditions. For the detailed relationships of these quantities it is referred to (Poblete et al., 2016).
The scalar function m is responsible for the stiffness increase after reversal loading: with the material parameter m R > 1. The scalar function y h = f(h,ε) introduces the influence of the Intergranular Strain h and is responsible for the reduction of the plastic strain rate upon cyclic loading. Hence, two threshold cases can be distinguished: • If y h = 0 the response of the model is rendered quasi-elastic and only negligible accumulation is obtained upon cyclic loading: • If y h = 1 the mechanical relation 6 reduces to the hypoplastic relation of von Wolffersdorff (1996) and the maximum plastic strain rate is rendered: The special feature of ISA (Fuentes and Triantafyllidis, 2015) compared to the Intergranular Strain concept (Niemunis and Herle, 1997) is the yield and the bounding surface within the Intergranular Strain space: whereby the material parameter R describes the maximum elastic strain amplitude and c is the hardening tensor describing the center of the yield surface. For further details of the constitutive relations it is referred to Poblete et al. (2016). The present implementation of the model uses a substepping scheme whereby small strain subincrements of approximately Δε ≈ 10 −6 are applied to achieve numerical convergence. Another yet noncommunicated modification applied in every implementation of ISA models for sands e.g. (Fuentes and Triantafyllidis, 2015;Poblete et al., 2016) is a so called "Phantom Elasticity" (PE). It is active during the complete calculation and provides an additional portion to the current stress stemming from Eq. (6): and the respective contribution to the numerical Jacobian. The isotropic elastic tensor C implies two additional material parameters: the Young's modulus E and the Poisson's ratio ν.
The modification for void ratios e exceeding the loosest possible void ratio e i (p) described in Section 5.1 was adopted. Furthermore, Appendix C addresses some additional differences between the hypoplastic version used in Poblete et al. (2016) compared to the version of von Wolffersdorff (1996).

Parameter calibration
Prior to the simulation of the vibratory pile driving tests, the constitutive parameters of the sand used in the experiments ("Karlsruhe sand") for the different constitutive models had to be determined 2 . The experiments used for the calibration of the material parameters consist of two oedometric compression tests (initially loose D r0 = 0.13 and dense D r0 = 0.96 samples) with loading, unloading and reloading (Fig. 4); three drained monotonic triaxial tests (initially loose D r0 = 0.26 and medium dense D r0 = 0.61 samples) with p 0 = 20 kPa and p 0 = 100 kPa and one undrained triaxial test with strain cycles of ε ampl = 0.05 (medium dense D r0 = 0.60 sample). The calibration procedure of each model will be summarised briefly in the following.

Hypoplasticity with Intergranular Strain
The parameters of the Hypo+IGS model have been calibrated as outlined in the following: • The eight parameters of the hypoplastic material model for "Karlsruhe Sand" have been calibrated based on loosely pluviated cones of dry sand (φ c ), index tests on maximum and minimum void ratio (e d0 , e c0 , e i0 ), oedometric compression tests (h s , n, β) and drained monotonic triaxial tests (α). In a first step, the parameters were calibrated by hand following the procedure proposed by Herle (1997).
• The hypoplastic parameters have been further optimised using the element test program Incremental Driver by A. Niemunis (Niemunis, 2008) to recalculate the laboratory tests. The parameters have been iteratively adjusted in order to receive a better fit of the material behaviour under monotonic loading conditions.
• The parameter R of Intergranular Strain (representing the elastic strain range) has been calibrated based on the normalised shear modulus (G/G 0 ) −γ curve presented in Rebstock (2011).
• The four parameters of Intergranular Strain m T , m R , β R and χ have been calibrated based on the undrained cyclic triaxial test. These parameters have been chosen with the aim to reproduce best the development of deviatoric stress amplitude and the accumulation of excess pore water pressure with increasing number of cycles measured in the cyclic triaxial test.
The parameters of the Hypo+IGS model used for the FE simulations are summarised in Table 2.

Sanisand
The 14 material parameters (except of the atmospheric pressure p a = 100 kPa and m = 0.05) of Sanisand have been calibrated using undrained as well as drained monotonic triaxial tests, oedometric compression tests and one undrained cyclic triaxial test. The parameters have been determined as follows: • The parameters of the critical state line e 0 , λ c and ξ have been initially determined by three undrained monotonic triaxial tests reaching the critical state and have been refined by the simulation of drained monotonic triaxial tests.
• M c and M e describe the critical state lines in the p −q plane for triaxial compression or extension, respectively. They have been determined by calculation using the Mohr-Coulomb relations To achieve better agreement with the results of the undrained cyclic triaxial test, M e was subsequently increased. Note that in the original model the parameter c = M c /M e is used as is done in Table 3.
• The parameters of the elastic material response (G 0 and ν) have been determined based on oedometric compression tests. In order to achieve comparable results with respect to the laboratory tests very low values were required for both parameters. This has also been noted in Wichtmann et al. (2019).
• The parameters h 0 , c h and n b have been calibrated based on the deviatoric stress -axial strain response of drained monotonic triaxial tests. Likewise, the parameters controlling the dilatancy (A 0 and n d )   2 To avoid confusion, it is again pointed out that different versions of the "Karlsruhe sand" exist in the literatur. In this context, a distinction has to be made in particular to the so-called "Karlsruhe fine sand' (d 50 = 0.14 mm, C U = 1.5)', which is much finer than the "Karlsruhe sand" used in this work. Parameters for "Karlsruhe fine sand" are presented e.g. in Wichtmann et al. (2019) and Liu et al. (2019).
have been determined based on the volumetric strain -axial strain curves of these tests.
• The parameters controlling the cyclic mobility effect (z max and c z ) have been calibrated based on the undrained triaxial test with strain cycles.
The final set of parameters is given in Table 3.

Hypoplasticity with intergranular strain anisotropy
For the calibration of the Hypo+ISA model parameters, the following steps were performed: • The eight parameters of the hypoplastic model have been calibrated analogously to the Hypo+IGS model (see Section 6.1).
• Similar to the Hypo+IGS model, the parameter R was initially chosen to fit the threshold strain amplitude (5⋅10 −5 ) which encloses the elastic locus of the Intergranular Strain from the (G/G 0 ) −γ curve presented in Rebstock (2011). However, to better reproduce the cyclic triaxial test, the parameter R was further reduced to R = 1⋅10 −5 .
• The minimum χ 0 and the maximum χ max Intergranular Strain exponent have been calibrated using the accumulated excess pore water pressure vs. the number of cycles curves from the undrained cyclic triaxial test.
• The accumulation rate factor C a controls the strain accumulation rate during the cycles. It has been iteratively adjusted in conjunction with the parameters m R and β R based on the undrained cyclic triaxial test.
• The two parameters of the numerically stabilising "Phantom Elasticity" have been adjusted to E = 40 kPa and ν = 0.45. These minimum values are indispensable in the present finite element calculation.
• All parameters have been further iteratively adjusted in order to receive the best possible fit of the material behaviour under  monotonic and cyclic loading with a unique parameter set. The bestfit parameters are summarised in Table 4.

Comparison based on laboratory tests
A respective comparison of experimental and numerical data is shown for the different constitutive models by means of recalculations of oedometric compression tests (Fig. 4), drained monotonic triaxial tests (Fig. 5) as well as an undrained cyclic triaxial tests (Fig. 6). It is noted that the parameter calibration was carried out with the aim of determining a single parameter set for each constitutive model that reproduces the different tests well on average. Thus, better agreement with the experiments could be achieved by determining a separate parameter set for each laboratory test. The simulations have been started at an initial vertical stress of σ v,0 = 1 kPa. In case of the Hypo+IGS and Hypo+ISA models, the Intergranular Strain has been assumed to be initially fully mobilised in the vertical direction. The parameters h s , n and β of the hypoplastic models (Hypo-+IGS and Hypo+ISA) have been calibrated to reproduce the first loading curves of the oedometric tests, thus the good prediction of these models at this phase of the test is expected. The oedometric stiffness upon unloading is also well reproduced by both hypoplastic models up to approximately σ v = 10 kPa. Afterwards (at small stresses) both hypoplastic models underestimate the stiffness and thus overestimate the vertical strain (deformation). This behaviour has been observed for all models (including the Sanisand model) and should be kept in mind for further constitutive model enhancements and developments.

Oedometric compression tests
All models show a poor performance during reloading for both the loose as well as the dense case since the reloading stiffness is strongly underestimated. For the Sanisand model, the comparison with the experimental results reveals that the oedometric stiffness is underestimated for all phases of the experiment (loading, unloading and reloading). The weaker performance of the Sanisand model under oedometric conditions, however, is known (see e.g. (Wichtmann et al., 2019)). An improvement of the performance of Sanisand for oedometric compression would be expected using the Sanisand version with closed yield-surface (Taiebat and Dafalias, 2008).

Drained monotonic triaxial tests
The simulations of three drained monotonic triaxial tests on airpluviated samples with different initial densities (D r0 = {0.26; 0.61}) and different initial mean effective stresses (p 0 = {20kPa; 100kPa}) are compared to the test results in Fig. 5.
On the top row of Fig. 5 the stress-strain curves are depicted. The initial stress was applied in the experiment by an isotropic increase of effective stress towards p 0 . For the hypoplastic models, a fully isotropic mobilisation of the initial Intergranular Strain has been assumed (h ii,0 = R/3). In case of Sanisand, the back-stress tensor α is initialised as the ratio of initial deviatoric stress tensor and mean effective stress. From the simulations, one may note that the peak strength is noticeably underestimated by all three models for the test with D r0 = 0.61 and p 0 = 100 kPa. For the sample with D r0 = 0.61 and p 0 = 20 kPa the Hypo+IGS and Hypo+ISA models slightly underestimate the maximum deviatoric stress q. Conversely, the Sanisand model predicts a slightly pronounced peak strength which slightly overestimates the one of the experiments. For the loose sample (D r0 = 0.26) all three constitutive models show an initially too stiff response but reproduce the residual strength fairly well. Generally, due to the almost identical constitutive equations (see the differences documented in Appendix C) governing the behaviour of the soil under monotonic loading, the Hypo+ISA model yields comparable simulation results with Hypo+IGS for the drained monotonic triaxial tests. It is worth noting that the disagreements of the q −ε 1 -curves could be reduced by further variation of the model parameters. However, given the loading conditions of the model tests, the emphasis of the calibration procedure was laid on the model response in a cyclic loading test and not on monotonic loading.
The volumetric strain behaviour ε v (ε 1 ) is displayed in the bottom row of diagrams in Fig. 5. For the loose sample (D r0 = 0.26) all three models show a good agreement with the experimental results. The ε v −ε 1 -curve of the dense sample at low mean effective stress (D r0 = 0.61, p 0 = 20 kPa) is well reproduced by the Sanisand model considering the possibility to control these curves almost independently (with the parameters n d and A 0 ) of the peak strength behaviour (p-q curves) once the parameters for the critical state line in the e −p space are fixed. For the experiment with D r0 = 0.61 and p 0 = 100 kPa, the Sanisand simulations shortcomings similar to those of the Hypoplastic ones. The inability of the hypoplastic models in reproducing the ε v (ε 1 ) behaviour of the dense samples is clearly visibly. This is because the parameter α controlling these curves was previously calibrated to reproduce the peak stress of the test.

Cyclic triaxial tests
The simulations as well as the experimental data of the undrained triaxial test with strain cycles of amplitude ε ampl = 0.05 are displayed in Fig. 6. The sample has been isotropically consolidated up to p 0 = 50 kPa reaching a medium initial density of D r0 = 0.6. This experiment was intentionally selected due to its low initial stress state, as is the case with the model tests of vibratory pile driving in the next section. Again, like in the monotonic tests, a fully isotropic mobilisation of the initial Intergranular Strain has been adopted (h ii,0 = R/3).
The effective stress paths depicted on the top row of Fig. 6 show a very good performance of the hypoplastic models. Even the deviatoric stress amplitude reached during the cycles is satisfactorily reproduced by the numerical simulations. On the other side, the Sanisand simulation underestimates severely the deviatoric stress in both triaxial compression and extension. The vanishing mean effective stress finally reached is well captured by all models. The decrease of the mean effective stress with the number of cycles results from the accumulation of the excess pore water pressure presented on the bottom row of Fig. 6. Obviously, all models reproduce the build-up of Δp w during the first 10 cycles as well as the final value of Δp w to the value of the initial mean effective stress. However, for N > 10, the rate of pore water pressure build-up is overestimated in all simulations.
A better performance of the simulation using Sanisand would be expected using the recently proposed extension by a memory surface as has been demonstrated in Liu et al. (2020,).

Comparison based on vibratory pile driving experiments
In the following, the capabilities of the chosen numerical approach in combination with the three constitutive models for the back-calculation of the model test described in Section 2 are investigated. The examination of the results includes comparisons of the pile displacements (Fig. 7), the pore pressure development at two different locations ( Fig. 11 and Fig. 12), the development of pile force vs. pile penetration at two different times of the experiment (Fig. 13 -Fig. 15) as well as contour plots of incremental vertical and horizontal displacements in the vicinity of the pile tip (Fig. 16).

Pile displacement
The development of the normalised pile displacement ũ y during the vibration is depicted in Fig. 7. The normalised pile displacement is defined as follows: with the calculated pile displacement u y and the pile diameter d Pile . The overall trend as well as the magnitude of pile displacement are reproduced well by all three constitutive models. Especially with regard to the complexity of the experiment, these simulation results are judged as satisfactory. The overall best agreement with the measurements is observed for the simulation with Hypo+IGS which delivered an almost perfect prediction of the pile displacement. However, some discrepancies are noticed. While the simulated and measured pile displacements during the first 3 s of pile driving are in nearly perfect agreement, the Sanisand model predicts a nearly constant pile penetration rate (for t > 3 s) and fails to capture the decrease in penetration rate observed in the experiment. This leads to an increasing deviation from the experiment resulting in an overestimation of final displacements of about 12%. Only small deviations from the experimental results are observed for the simulation with the Hypo+ISA model. Compared to the other two simulations, the Hypo+ISA model predicts the strongest decrease in pile penetration rate with ongoing vibratory driving.
The displacement-time histories of the experiment and the simulations are illustrated in detail for the time frames S1, S2 and S3 in Fig. 8. The comparison shows that the amplitudes in the simulations during time frame S1 are too large compared to the measured ones in the experiment (although the overall displacement trend is reproduced well, see also Fig. 7). In addition, a small phase shift (indicated by the vertical lines in Fig. 8) between the simulations and the experiment is observed. In time frame S2, both the displacement trend as well as the displacement amplitudes in the simulations and the experiment are in good agreement. Similar observations are made for time frame S3. The initially small phase shift increases with ongoing pile penetration and becomes most apparent in time frame S3.

Pile acceleration
A comparison of the measured and calculated (vertical) pile accelerations is provided in Fig. 9. In the experiment, the acceleration magnitude during the downward movement of the pile (negative values) exceeds the one of the upward movement (positive values). After the first second of pile driving, the double acceleration amplitude is approximately 50 m/s 2 with a tendency to increase towards the end of the pile driving process. Contrary to the experiment, the differences in acceleration magnitude during the downward and upward movement of the pile are absent in the simulations. Only for the simulation with the Hypo+IGS model a slight increase in magnitude during the downward movement of the pile is observed with increasing pile penetration, however less pronounced than in the experiment. Compared to the experiment, the simulations with the Sanisand and Hypo+IGS model slightly overestimate the acceleration amplitude. The simulation with the Hypo+ISA model shows the best agreement with the measured acceleration magnitudes.
The corresponding response (Fourrier) spectra of the accelation-time histories are given in Fig. 10. The decisive frequency of 25 Hz (excitation frequency of the vibrator) is clearly visible in the experiment as well as in the simulations. However, as for the acceleration-time histories, the magnitude is overestimated in the simulations. For higher frequency Fig. 8. Comparison of normalised pile penetration according to Eq. (11) over vibration time for the time frames S1, S2 and S3 marked in Fig. 7.   Fig. 9. Comparison of pile acceleration over vibration time for the experiment (black) and the simulations (coloured). components, this relationship reverses. Due to the numerical damping of higher frequency components in the simulations (introduced by the HHT time integration scheme, see Section 3), the high frequency components are damped too much (compared to the experiment).

Excess pore water pressure
Changes in pore water pressure were measured at two pore pressure transducers PPT A and PPT B (see Fig. 1). The comparisons of the results of the simulations with the measurements are provided in Fig. 11 and 12 for PPT A and PPT B, respectively. Therein, Δp w denotes the change in pore water pressure and p w 0 is the initial (hydrostatic) pore water pressure. Referred to the driving process in the experiment, PPT A is passed by the pile tip after the first seconds of vibration whereas PPT B is reached after approximately 6 s. Comparing the measured development of pore water pressure for PPT A to the calculated ones (Fig. 11), one may conclude that, although a reasonable agreement between the experiment and simulations is observed, all three models show different shortcomings. For the simulation with the Hypo+IGS model a good agreement of the minimum values of pore water pressure (Δp w /p w 0 < 0, i.e. decrease with respect to the hydrostatic values) during vibration is observed. However, the maximum values of pore water pressure (and thus the tendency towards contraction under cyclic loading) during the first two seconds of pile driving are strongly underestimated in the simulation. The opposite is observed for the simulation with the Hypo+ISA model: while the increase in pore water pressure is captured satisfactory, the Hypo+ISA model fails to reproduce the decrease in pore water pressure and thus does not adequately capture the tendency towards dilatation under cyclic loading. Contrary to the simulations with the hypoplastic models, the simulation with the Sanisand model is able to capture both, positive as well as negative pore pressure increments. However, the Sanisand simulation drastically overestimates the change in pore water pressure during the first two seconds. When the pile tip has passed the transducer (approximately after 2 s) the amplitude of incremental pore water pressure is underestimated in all simulations.
Similar observations are made for the comparison of measured and calculated excess pore water pressure development at transducer PPT B. In the experiment a slow increase in pore pressure amplitude is observed during the first four seconds of vibratory driving. Then, as the pile tip enters the vicinity of the transducer, a strong increase of the amplitude of pore water pressure takes place. This slow increase of Δp w during the first four seconds is captured well by the simulation with the Hypo+IGS model. The simulations with the Sanisand and Hypo+ISA models, however, predict an almost constant amplitude and an accumulation of excess pore water pressure during this phase of the simulation. The strong increase of the amplitude of excess pore water pressure starting at t ≈ 4 s is apparent in all three simulations. However, similar to the observations for PPT A, the magnitude of Δp w is overestimated in the simulation with the Sanisand model and underestimated in the simulations with the Hypo+IGS and Hypo+ISA models, respectively. The decrease in pore water pressure amplitude after about 5.5 s in the  Fig. 11. Development of pore water pressure recorded at PPT A (see Fig. 1) for the experiment (black) and the simulations (coloured). simulation with the Sanisand model is a consequence of the overestimated pile penetration rate (see Fig. 7): the maximum pore water pressure is reached when the pile tip reaches PPT B. As soon as the pile tip has passed the pore pressure transducer, a decrease of pore water pressure takes place (compare Fig. 11). It should be noted that in the simulation with the Sanisand model oscillations in pore water pressure are observed even during step 4 of the simulation when the pile is at rest. These oscillations are absent in the experiment and the other simulations.
In the simulations with the Hypo+IGS model and the Sanisand model some numerical instabilities by means of sharp increases (or decreases) of excess pore water pressure are observed. For the simulation with the Hypo+IGS model, these instabilities are apparent only during the first second of the simulation.
In the authors' opinion, the too large pore water pressure amplitude and the overestimated penetration rate of the pile in the simulation with the Sanisand model can be traced back to shortcomings observed in the back-calculation of the laboratory tests: the predicted initial loading stiffness in the oedometric compression tests was too low and a too fast relaxation of effective stress in the undrained cyclic triaxial test was observed. It is all the more surprising that noticeable differences were found in the simulations with Hypo+IGS and Hypo+ISA, although very similar results were obtained in the recalculation of the laboratory tests.

Pile force vs. pile displacement
To further evaluate the performance of the different constitutive models for the prediction of the pile penetration behaviour, the measured and calculated pile force as a function of the normalised pile displacement ũ y is compared in Fig. 13. The normalised pile force reads: Therein, the pile inertia forces m Pile ⋅a Pile y are subtracted from the pile head forces F Pile (measured between the unbalances and the pile head) during the vibratory pile driving. The resulting force is then related to the static pile force F Pile stat resulting from the dead weight of the pile, the oscillator and the load cell to achieve comparability between the simulation (full model) and the experiment (nearly half model). Accounting for the model dimensions, the pile mass is m Pile = 0.627 kg in case of the experiment (nearly half model) and m Pile = 1.085 kg for the simulation (full model). The pile force F Pile is evaluated at a load cell Fig. 13. Normalised pile force F Pile as function of the normalised pile displacement ũ y .

Fig. 14.
Normalised pile force F Pile as function of the normalised pile displacement ũ y for time frame S1 (see Fig. 13).

Fig. 15.
Normalised pile force F Pile as function of the normalised pile displacement ũ y for time frame S2 (see Fig. 13).
installed between the pile head and the vibrator. Both the measured and the calculated pile head force F Pile contain the skin friction as well as the tip resistance. In fact, the development of the pile force in the simulations reflects the mobilised soil resistance during pile driving.
In the experiment, an increase in normalised pile force F Pile (both positive and negative) with increasing pile penetration ũ y is observed. While the pile force at the beginning of the pile driving is in the range −0.6⩽F Pile ⩽1.5, this range grows with ongoing penetration depths until it reaches final values of −1.0⩽F Pile ⩽3.5. This steady increase in pile force is not captured well in the simulations. While the simulation with the Hypo+IGS model shows an increase in pile force (both positive as well as negative), this increase is not as steady and pronounced as in the experiment. Contrary to the experiment, the simulation with the Sanisand model predicts a nearly constant pile force with increasing pile penetration which fits to the predicted constant pile penetration rate observed in Fig. 7. Both, the simulations with the Hypo+IGS and the Sanisand model start from pile forces similar to the experiment at the beginning of the vibratory pile driving. However, with increasing pile penetration, the magnitude of pile force is underestimated in both simulations. The simulation with the Hypo+ISA model shows the least agreement with the experiment.
A more detailed comparison of the measured and calculated pile forces for two time frames S1 and S2 are provided in Fig. 14 and Fig. 15, respectively. Overall, all simulations reproduce the measured development of pile force qualitatively: During the penetration phase (Phase 1-3) a significant increase in pile force is observed. During the upwards movement of the pile (Phase 3-5), the normalised pile force decreases rapidly and eventually reaches a small negative value, which remains approximately constant during the rest of the upwards movement. This negative force mainly results from friction between the pile shaft and the soil. This is more pronounced during time frame S2 (Fig. 15) than during time frame S1 (Fig. 14).
For time frame S1 (Fig. 14), the development of pile force is most accurately captured by the simulation with the Hypo+IGS model. The simulation with the Sanisand model captures the qualitative development of pile force reasonably well, but overestimates the penetration rate. Both, the simulations with the Hypo+IGS and the Sanisand models underestimate the maximum pile force slightly, whereas the simulation with the Hypo+ISA predicts noticeably too small maximum pile forces. The sharp decrease of pile force during the first part of the upward movement of the pile (Phase 3-4) is underestimated in all simulations. In addition, the constant force during the upward motion of the pile is only partly captured by the simulations. Contrary to the experiment, an irregular development of pile force is observed in the simulation with the Hypo+ISA model.
During time frame S2 (Fig. 15), the magnitude as well as the development of pile force are again most accurately captured by the simulations with the Hypo+IGS and the Sanisand models. In the experiment, an increase in maximum pile force of about 30% and a decrease of permanent pile displacement during one cycle compared to time frame S1 are observed. Both characteristics are qualitatively captured by these simulations. However, the increase in pile force is only of about 20% for the simulation with Hypo+IGS model and smaller than 10% for the simulation with the Sanisand model. As for the simulation with the Hypo+ISA model, the irregularities observed during time frame S1 are absent during time frame S2. However, the increase in pile force with  increasing penetration is not captured by this simulation. Moreover, the simulation with the Hypo+ISA model even predicts a decrease of pile force during the second half of the downward movement of the pile (Phase 2-3).

Incremental displacement fields
The analysis of the model tests also comprises spatial distributions of displacement and strain in the vicinity of the pile tip using digital image correlation (DIC). A detailed documentation of how the results were obtained and evaluated is given in Vogelsang et al. (2016) and Vogelsang (2017). For the present comparisons the incremental horizontal and vertical displacement during the downwards movement of the pile within one cycle in time frame S3 are analysed. The evaluation of incremental fields is helpful to assess the performance of the constitutive models since it allows a local comparison with the experimental results opposite to the previous global comparisons (pile displacement and driving force). Fig. 16 a,b,c display the change of horizontal displacement in dependence of the spatial location during the pile movement Δu Pile y illustrated in Fig. 16 a (phase 1-2 according to Fig. 14). Each plot displays the spatial distribution observed in the experiment on the left-hand side and the results of the simulations using the different constitutive models on the right-hand side. The location of the pile tip is indicated by the schematic sketch (hatched). The incremental horizontal displacement Δu x is divided by the increment of vertical pile displacement Δu Pile y during the considered phase of the pile movement. Note that Δu Pile y is different for each simulation and is given at the top of each comparison. The simulation using the Hypo+IGS model displayed on the right-hand side of Fig. 16a  ⃒ are substantially smaller than the ones obtained from the measurements. The simulation using the Hypo+ISA model again shows similar results as the simulation using the Hypo+IGS model. In the vicinity of the pile, however, the Hypo+ISA model represents the experimental results worse.
The evaluation of the incremental fields shows that all three constitutive models are capable to adequately represent the local soil behaviour. The Hypo+IGS and Hypo+ISA model predict similar incremental deformations whereas the simulation using Sanisand shows a greater spreading of vertical incremental displacement.

Summary and conclusions
Laboratory tests (oedometer, monotonic and cyclic triaxial tests) and model tests of vibratory pile driving in saturated "Karlsruhe Sand" were simulated using three different constitutive models. By comparing the simulation results with the experimental results the performance of the numerical implementation and of the soil models was evaluated. The parameter calibration was carried out with the aim of determining a single parameter set for each constitutive model that reproduces the different tests well on average. The capabilities of the constitutive models to reproduce the laboratory tests can be summarised as follows: • The oedometer tests were best reproduced by the Hypo+IGS and the Hypo+ISA models. Both models showed a good agreement with the experiment on loose and dense samples during the first loading and unloading cycle. However, the reloading stiffness was underestimated. The Sanisand model was not able to predict oedometric tests properly with the chosen parameter set.
• The drained monotonic triaxial tests were best reproduced by the Sanisand model. Contrary to the hypoplastic constitutive models, both the q −ε 1 -and ε v −ε 1 -curves were captured satisfactory. The Hypo+IGS and Hypo+ISA models showed limitations in reproducing the ε v −ε 1 -curves.
• The undrained cyclic triaxial test was best fit by the Hypo+IGS and the Hypo+ISA models, which out-performed the Sanisand simulations in capturing the decrease of the deviatoric stress amplitudes with increasing effective stress relaxation. All three models successfully predicted the accumulation of excess pore water pressure. However, compared to the experiment the build-up of pore water pressure occurred too fast.
The conclusions drawn from the simulation of the element tests are in agreement with a similar investigation on another sand reported in Wichtmann et al. (2019).
The performance of the three constitutive models was further investigated by the back calculation of vibratory pile driving model tests in saturated sand. It was shown that a parameter calibration based on only a few, but carefully chosen, laboratory tests can lead to satisfactory simulation results. The main findings are: • In combination with the numerical implementation, all three constitutive models were able to reproduce satisfactorily the pile displacement (penetration) during the vibratory driving process. However, the simulation with the Sanisand model underestimated the decrease of the pile penetration rate and thus overestimated the final pile penetration depths by approximately 12%.
• The development of excess pore water pressure during the vibratory pile driving was captured acceptable by all three constitutive models. However, all models showed some weaknesses in accurately predicting the magnitude of excess pore water pressure amplitude.
• The development of pile force with increasing pile penetration could not be accurately reproduced by the simulations. However, an acceptable agreement was achieved by the simulation with the Hypo+IGS model, followed by the simulation with the Sanisand model. The Hypo+ISA model could neither reproduce the development of pile force qualitatively nor quantitatively. The reason for the good performance in predicting the pile displacement but simultaneously failing in capturing the pile force using the Hypo+ISA model could not yet been fully understood and is subject of further investigations.
• In the experiments, incremental displacement fields have been monitored using DIC. The best agreements for these fields were observed for the simulations with the Hypo+IGS and Hypo+ISA models. A slightly worse, but still quite acceptable agreement was found for the simulation with the Sanisand model.
The shortcomings of the constitutive models observed in the simulations of the element tests were partly also apparent in the simulation of the vibratory pile driving model tests. For instance, the Sanisand model predicted a too low initial loading stiffness in the oedometric compression tests and a too fast relaxation of effective stress in the undrained cyclic triaxial test. Similar shortcomings were present in the simulation of the model test as the penetration rate of the pile, the pore-water pressure amplitudes and the incremental vertical displacement were found to be too large. Some of these issues could be solved by increasing the elastic material constants (G 0 and ν). However, severe numerical instabilities for the simulation of the vibratory pile driving model tests were observed using a set of parameters with adjusted elastic properties. Future work will concentrate on resolving these issues and further investigating the influence of the set of parameters on the performance of the constitutive models in BVPs.
Overall, the present study demonstrates that existing advanced constitutive models can be successfully applied to the numerical investigation of complex geotechnical BVPs. The investigated BVP places high demands on the numerical methods and the constitutive models as it faces high-frequency cyclic loading of saturated soils at small stress states. Especially the cyclic shearing in the contact zone of the pile and the soil require numerically stable and robust implementations. All three constitutive models were able to meet these requirements. Especially the good agreement of the predicted pore water pressures and incremental displacement fields enables the investigation of further questions, such as the influence of the pile installation on neighbouring piles.

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.
The objective stress rate σis calculated from σ = M :ε, where ε is the strain rate. The stiffness tensor M is defined by Niemunis and Herle (1997) where β R is another material parameter and I is the fourth-order identity tensor. The stiffness tensors L and N are given by the following equations (von Wolffersdorff, 1996): Therein σ = σ trσ and σ * =σ − 1 3 1 are used. For information on the factor f e i it is referred to Section 5. φ c , h s , n, e i0 , e d0 , e c0 , α and β are parameters and e is the actual void ratio. The pressure-dependent void ratios e i , e c and e d in Eq. (A.8), describing the loosest, the critical and the densest state, are calculated using the following relation (Bauer, 1996)

Appendix B. Equations of the Sanisand model
In the following, the implemented equations used for the simulations with the Sanisand model are given. They are identical to those presented in Dafalias and Manzari (2004). The objective stress rate is calculated using: The elasto-plastic stiffness E ep is: The change of the fabric-dilatancy tensor z with plastic volumetric strain increment is Lastly, the change of the back-stress tensor α is The initial back-stress tensor α ini is updated to α in case of a load reversal which occurs in case of (α −α ini ) : n < 0. B and C in Eq. (B.6) are: is the convective coordinate alongside the edge of the surface i = 2, N J is the shape function of node J of the surface i = 2, x J and x (1) I are the global coordinates of node J of surface i = 2 and node I of surface i = 1, respectively. x (2) ,ξ (ξ (2) ) is the derivative of the global coordinates of surface i = 2 with respect to ξ evaluated at ξ (2) . As indicated in Fig. D.17, ξ is the local coordinate of the surface point at which the minimum distance to a node of the paired surface is found. Eq. (D.1) is evaluated for the corner and the middle nodes of both surfaces. Based on ξ, the minimum distance between the node I, to which the projection is performed, and surface 2 is found.
Eq. (D.1) is solved iteratively using Newton's method. Once converged, the normal gap between the surfaces is calculated using ⋅x (2) ,ξ .

(D.2)
The normal contact traction t N = εg N is calculated using the penalty factor ε. Note that the normal contact traction is the total traction in case of the u-p formulation. A Coulomb friction model with non-associated flow rule is used (no dilatancy is accounted for). The effective normal traction t ′ N , calculated by subtracting the pore water pressure from the total normal traction at every node (or by interpolation in case of the middle nodes) is used to calculate the (in 2D) scalar tangential traction t T . The tangential traction is calculated incrementally using Δt T = E T Δg T . E T is the tangential shear stiffness which can be obtained by interface shear tests. The tangential traction is restricted by Eq. (D.3), where δ is the wall friction angle: For the simulation of the vibratory pile driving a tangential stiffness of E T = 3000 kPa and a wall friction angle of δ ≈ 1/3φ c = 11 • are used. The first value has been obtained from measurements in interface tests reported in Vogelsang (2017), while the latter is an estimate considering the smooth surface of the pile.