Accepted Manuscript

Stone column settlement performance in structured anisotropic clays: the influence of creep

Brian G. Sexton, Bryan A. McCabe, Minna Karstunen, Nallathamby Sivasithamparam

PII: S1674-7755(16)30046-4

DOI: 10.1016/j.jrmge.2016.05.004

Reference: JRMGE 260

To appear in: Journal of Rock Mechanics and Geotechnical Engineering

Received Date: 23 February 2016 Revised Date: 14 April 2016 Accepted Date: 19 May 2016

Please cite this article as: Sexton BG, McCabe BA, Karstunen M, Sivasithamparam N, Stone column settlement performance in structured anisotropic clays: the influence of creep, Journal of Rock Mechanics and Geotechnical Engineering (2016), doi: 10.1016/j.jrmge.2016.05.004.

Rock Mechanics and

Geotechnical

Engineering

—- —

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A rrCDTCn N/fA MI TC nT> TÜT

Journal of Rock Mechanics and Geotechnical Engineering

Journal online: www.rockgcotteh.org

Stone column settlement performance in structured anisotropic clays: the influence of creep

Brian G. Sexton1, Bryan A. McCabe1, , Minna Karstunen2, Nallathamby Sivasithamparam3

1. College of Engineering and Informatics, National University of Ireland, Galway, Ireland

2. Department of Civil and Environmental Engineering, Chalmers University of Technology, Gothenburg, Sweden

3. Computational Geomechanics Division, Norwegian Geotechnical Institute, Oslo, Norway Received 23 February 2016; received in revised form 14 April 2016; accepted 19 May 2016

Abstract: The recently developed elasto-viscoplastic Creep-SCLAY1S model has been used in conjunction with PLAXIS 2D to investigate the effectiveness of vibro-replacement in a creep-prone clay. The Creep-SCLAY1S model accounts for anisotropy, bonding, and destructuration, and uses the concept of a constant rate of viscoplastic multiplier to calculate creep strain rate. A comparison of settlement improvement factors with and without creep indicates that 'total' settlement improvement factors (primary plus creep) are lower than their 'primary' counterparts (primary settlement only). The lowest settlement improvement factors arise for analyses incorporating the effect of bonding and destructuration. Examination of the variations of vertical stress with time and depth has indicated that vertical stress is transferred from the soil to the column as the soil creeps. This results in additional column yielding. In addition, the radial and hoop stresses in the soil are lower for the 'creep' case. The reduced radial stresses lead to additional column bulging and hence more settlement, whereas the hoop stress reductions appear to be a secondary effect, caused by additional plastic deformation for the 'creep' case.

Keywords: stone columns; creep; anisotropy; destructuration; finite element (FE) method

1. Introduction

Vibro-replacement has traditionally been considered to be an effective means of improving the bearing capacity and settlement characteristics of mixed fills and weak soils. Vibro-replacement solutions can be more cost-effective, quicker to implement and less CO2-intensive than piling alternatives in certain circumstances. The technique is becoming increasingly popular for the treatment of soft fine-grained deposits; some of these soils can comprise a significant organic content, in which case creep settlements can contribute a significant proportion of the total settlement.

The settlement reduction potential of vibro stone columns is typically quantified using a settlement improvement factor (n):

n = ¿untreated/¿treated (1)

where 5untreated and ¿Med are the final settlements of untreated (i.e. no columns) and treated ground, respectively.

The majority of analytical settlement design methods pertain to primary settlement only (e.g. Priebe, 1995; Castro and Sagaseta, 2009; Pulko et al., 2011), and so the same n value tends to be applied to both primary and creep settlements in routine designs. In addition, the majority of n values measured in the field (McCabe et al., 2009) tend to be 'lumped', with no distinction between initial compression, primary consolidation settlement, and creep. The length of time required to measure 'pure' long-term creep settlements in soft low permeability soils serves as the main impediment in the latter case.

While the majority of laboratory-scale testing carried out to date has been informative (e.g. Black et al., 2011), it tends to be limited by scale effects and a difficulty in replicating realistic boundary conditions. Additionally, Castro et al. (2013) have noted that the reconstituted soils

* Corresponding author. Tel.: +353 91492021. E-mail address: bryan.mccabe@nuigalway.ie

used in laboratory testing are not fully representative of natural clay behaviour, while Mesri (1973) has pointed out that the rate of creep is lower in reconstituted soils. There are also difficulties associated with extrapolating long-term performance in the field from short-term laboratory tests (Mitchell and Kelly, 2013).

Moorhead (2013) carried out a series of laboratory tests in one- (1D) and three-dimensional (3D) loading chambers, and is, to the authors' knowledge, the only researcher to date to investigate the creep settlement reduction potential of vibro stone columns in the laboratory. The tests were carried out on reconstituted samples of kaolin and Belfast sleech and examined stone column behaviour for both a rigid raft and an isolated loading scenario. Although the laboratory data showed a significant amount of scatter, it was tentatively concluded that columns effectively reduced primary settlement at low bearing pressures but were ineffective at high pressures, and had only a minor influence on reducing initial compression and creep settlements. The findings need to be treated with caution because the initial conditions in the untreated and treated soil beds were different in some cases.

In this paper, a series of axisymmetric analyses carried out using the PLAXIS 2D finite element (FE) program (Brinkgreve et al., 2011) is reported with a view to assess the settlement reduction potential of vibro-replacement in a structured anisotropic creep-prone clay. The Bothkennar soft clay test site in Scotland, comprising an overconsolidated crust overlying two layers of lightly overconsolidated Carse clay, has been used as a representative soil profile for the numerical analyses.

In previous studies, Sexton and McCabe (2013, 2015) carried out some preliminary numerical work using a simplified single-layer version of the multi-layer Bothkennar profile to gauge the influence of creep on settlement improvement factors. Subsequently, Sexton and McCabe (2016) extended this work to a multi-layer scenario. The commercially available soft soil creep (SSC) model was used to represent the host clay behaviour in these studies. The latter study in particular provided valuable insight into the likely behaviour of stone columns in creep-prone soils and

formed an important frame of reference for one using a more advanced constitutive model, such as the Creep-SCLAYIS model, which is used to model the soft clay in this paper. The Creep-SCLAYIS model (Sivasithamparam et al., 2013, 2015) incorporates anisotropy, bonding, and destructuration, each of which can be 'switched off individually or in various combinations by adjusting the input parameters. The model is not yet commercially available, and therefore a user-defined model implementation into the PLAXIS FE code was used. The hardening soil (HS) model (Schanz et al., 1999) is used to represent the granular column material. For simplicity, any installation effects have been ignored.

2. Modelling creep using advanced constitutive models

2.1. Model classification

Constitutive models for describing the time-dependent behaviour of soft soils can be classified as either empirical models, rheological models, or general stress-strain-time models, each of which have been reviewed in detail by Liingaard et al. (2004). Empirical models are generally obtained by fitting mathematical/constitutive expressions to experimental data whereas rheological models tend to be used to gain a conceptual understanding of time effects in soil.

General stress-strain-time models are capable of describing the rate-dependent behaviour of soils under a variety of different loading conditions. These models tend to be formulated incrementally and so are suitable for implementation within the FE method. The majority of elasto-viscoplastic general stress-strain-time models are based on overstress theory (e.g. Perzyna, 1963, 1966), either 'conventional overstress' or 'extended overstress'. 'Extended overstress models' are preferable to 'conventional overstress models' (Yin et al., 2010).

2.2. 3D elasto-viscoplastic models

2.2.1. Isotropic models

The commercially available isotropic SSC model (Vermeer et al., 1998; Vermeer and Neher, 1999) can be classified as either an 'extended overstress model' or a 'creep model', as defined by Yin et al. (2010). 'Creep models' use the creep coefficient, Ca, or its isotropic equivalent, n', as the soil viscosity input parameter. The isotropic elasto-viscoplastic model developed by Yin et al. (2002) is also denoted a 'creep model'.

2.2.2. Anisotropic models

Anisotropic elasto-viscoplastic soil models have been developed as extensions to the EVP and SSC models, e.g. the anisotropic EVP model (Zhou et al., 2005) and the anisotropic creep model (ACM, Leoni et al., 2008), respectively. These 'creep models' assume that the viscoplastic volumetric strain rate is independent of the stress state, and consequently they predict unrealistic strain-softening behaviour for undrained triaxial tests on isotropically consolidated samples (Yin et al., 2010; Sivasithamparam et al., 2013, 2015). To overcome this deficiency, Yin et al. (2010) proposed a new anisotropic elasto-viscoplastic soil model in which the volumetric strain rate is dependent on the stress ratio, n = q/p' (where p' and q denote the mean effective stress and deviatoric stress, respectively).

2.2.3. Anisotropic models with bonding and destructuration

The Creep-SCLAY1S model (Sivasithamparam et al., 2015) is an anisotropic soil model that also takes account of bonding and destructuration. 'Destructuration' refers to the progressive breakdown/degradation of bonds during straining (Leroueil et al., 1979) and is accommodated using the concept of an intrinsic yield surface proposed by Gens and Nova (1993). Other models which fit in this category are the AniCreep model developed by Yin et al. (2011) as an extension to the model developed by Yin et al. (2010), the EVP-

SCLAY1S model (Yin and Karstunen, 2008) which is categorised as a 'conventional overstress model', and the non-associated structured anisotropic creep (n-SAC) model developed by Grimstad et al. (2010).

The key feature of the n-SAC model is that the time resistance concept is introduced to the viscoplastic multiplier rather than to the viscoplastic volumetric strain. Models which assume constant contours of volumetric creep strain (e.g. SSC or ACM) yield unrealistic creep strains for almost all stress paths (Olsson, 2013). In keeping with the n-SAC model, the Creep-SCLAY1S model also uses the concept of a constant rate of viscoplastic multiplier to calculate creep strain rate, but the commonly used semi-logarithmic creep coefficient, Ca, is used as the viscosity input parameter, making the model attractive from a practical modelling perspective.

3. Formulation of the Creep-SCLAY1S model

The formulation of the Creep-SCLAY1S model (in triaxial stress space) was described by Sivasithamparam et al. (2013) and it was extended to 3D for FE analyses by Sivasithamparam et al. (2015). In the version used in this paper, destructuration has also been incorporated. The anisotropy and destructuration components of the model were formulated using the constitutive equations from the rate independent S-CLAY1 (Wheeler, 1997; Wheeler et al., 2003) and S-CLAY1S (Koskinen et al., 2002; Karstunen et al., 2005) models, which accounted for anisotropy and both anisotropy and destructuration, respectively.

For simplicity, the Creep-SCLAY1S model is briefly explained here in triaxial stress space. For the extension to 3D, readers can refer to Sivasithamparam et al. (2015). The total strain rate (e) is composed of an elastic component ( ee ) and a viscoplastic (creep) component ( ec), i.e. e = £' + ec (2)

The rotational hardening law, describing the changing inclination of the yield surface due to creep strains, takes the form shown as follows:

{j - a) Kc) + wd (j -aj N

where w and wd are two additional soil constants; de£ and ded are the increments of creep volumetric and deviatoric strains, respectively; and () denote Macaulay brackets. The constant wd controls the relative effectiveness of the deviatoric and volumetric creep strains in determining the overall target value for a, while w controls the absolute rate at which a approaches the target value, where a is the angle of inclination of the yield surface (Fig. 1).

Fig. 1. Yield surfaces of the Creep-SCLAYl S model in triaxial stress space.

The destructuration law describing the breakdown of bonding caused by creep strains takes the following form:

dx=-xx{ M +XdN ) (4)

where £ and are two additional soil constants controlling the absolute rate of destructuration and the relative effectiveness of volumetric and deviatoric creep strains, respectively, in destroying the bonding.

The initial amount of bonding (x0) relates the size of the natural yield surface ( pp0 ) to the size of the intrinsic yield surface ( p'p0i), i.e. P'po = (1 + Xo) Ppoi (5)

The yield surface (normal consolidation surface, NCS) evolves with the creep volumetric strains according to Eq. (6), with the equivalent mean stress measure ( p'eq , Eq. (7)) defining the intersection of the current stress surface (CSS) with the p' axis (Fig. 1), i.e.

Pp = Ppo exP

Peq = P

l -к (q -ap')2

[M 2(ff) -a2]p'

where М(в) is the stress ratio at critical state, which is a function of Lode angle (в) to incorporate a smooth failure surface (see Sivasithamparam et al., 2015); and A" and к* are the modified compression and swelling indices, respectively, and hence related to the 1D compression and swelling indices (Cc and Cs).

In contrast to the SSC model and the ACM, which calculate the volumetric creep strain rate according to Eq. (8), the Creep-SCLAY1S model uses the rate of viscoplastic multiplier (Л) (Eq. (9)):

Л = -

" m \в) -ol

M г(в)-h

where the additional term in Eq. (9) was added to ensure that under oedometric conditions, the resulting volumetric creep strain corresponds to Eq. (8); a nc denotes the inclination of the yield ellipse in the normally consolidated condition; and t is a reference time, which is usually 1 day if n' is calculated using an incremental load oedometer test with a loading duration of 1 day.

The determination of the additional model parameters required for the Creep-SCLAY1S model is straight-forward. The anisotropy parameters,

Table 1. Bothkennar

a0 (initial yield surface inclination) and wd can be calculated from the critical state friction angle ( ф ) and the coefficient of lateral earth pressure in the normally consolidated condition (K0nc). The value of w should be calculated/optimised by simulating undrained triaxial extension tests, or in their absence, simply estimated based on compressibility (Zentar et al., 2002).

The initial amount of bonding, x0, can be calculated based on the sensitivity, St (see Eq. (10)), and the other destructuration parameters and £) can be calibrated using the optimisation procedure described in Koskinen et al. (2002). The intrinsic compression and creep indices, ll and mm (measured from oedometer tests on reconstituted samples), should be used as opposed to A* and ц' (measured from oedometer tests on natural samples) when modelling destructuration using the FE method. X = S, -1 (10)

4. Soil profile

4.1. Soil parameters

The Bothkennar soft clay test site in Scotland was purchased by the UK Science and Engineering Research Council (SERC) in 1987 as a national soft clay test bed. The silty clay at Bothkennar is highly structured with an organic content of 3%-5%, depending on the 'facies' type (bedded, laminated, mottled and weathered) (e.g. Paul et al., 1992), and a bulk unit weight of у = 16.5 kN/m3 (e.g. Nash et al., 1992a). The multi-layer soil profile (Table 1) adopted in this study is based on the HS model soil profile used by Killeen and McCabe (2014). The authors obtained the parameters from ICE (1992) and validated their profile against a field load test on an unreinforced rigid pad footing described by Jardine et al. (1995). However, creep was not considered in their study and so the additional creep parameters in Table 1 were calculated based on a Ca/Cc (= ц /A ) ratio of 0.04 (e.g. Nash et al., 1992b). The additional anisotropy and destructuration parameters for Bothkennar clay quoted in Table 2 have previously been calibrated by Karstunen et al. (2013).

The initial stress state for the FE model has been generated using a pre-overburden pressure (POP) of 15 kPa for the upper layers and an overconsolidation ratio (OCR) of 1.5 for the lower Carse clay, with S0 and denoting the initial effective stress and 1D preconsolidation stress, respectively. The in situ at-rest earth-pressure coefficients (K0) are based on a series of spade cell, self-boring pressuremeter (SBPM), and Marchetti dilatometer tests carried out by Nash et al. (1992a).

material parameters.

Type Depth (m) у (kN/m3) OCR POP (kPa) Ko eo A* A,* к И И* c' (kPa) V (°) V kx (m/d) ky (m/d)

Crust 0-1.5 18 - 15 1.5 1 0.015 0.006 0.002 0.0006 0.0002 3 0 0.2 1 x 10-4 6.9 x 10-5

Upper Carse clay 1.5-2.5 16.5 - 15 1 1.2 0.049 0.018 0.006 0.002 0.0007 1 0 0.2 1 x 10-4 6.9 x 10-5

Lower Carse clay 2.5-14.5 16.5 1.5 - 0.75 2 0.162 0.06 0.023 0.0065 0.0024 1 0 0.2 1 x 10-4 6.9 x 10-5

Note: OCR = S / SO ; POP = sP - s0

Table 2. Additional material parameters.

Mc Me a ffld w X0 id i

1.375 -1 0.5267 0.9281 50 8 0.2 9

The slopes of the critical state lines (CSLs) in compression (Mc) and extension (Me) have been selected based on a series of triaxial stress paths tests on reconstituted clay carried out by Allman and Atkinson (1992). The Mc corresponds to a critical state friction angle of 34°; this high friction angle is attributable to both a high organic content and an abundance of silt-sized grains. Nominal cohesion values (c' = 1 kPa) have been used for numerical stability and a dilatancy angle of 0° was

adopted as representative of a lightly overconsolidated clay. The horizontal and vertical permeabilities (kx and ky) were measured by Leroueil et al. (1992) using both in situ (e.g. pushed-in-place piezometers, self-boring permeameters, BAT system) and laboratory (e.g. oedometer cells, triaxial cells, radial flow cells) methods.

The adopted soil parameters derived predominantly from the results of oedometer tests have been validated by using the PLAXIS 'Soil Test' facility to simulate the undrained triaxial compression tests reported by Atkinson et al. (1992) (see Sexton, 2014). 4.2. Scenarios

Three separate scenarios have been considered for the Creep-SCLAY1S model analyses described in Section 6:

(1) Anisotropy and destructuration have been 'switched off by setting the relevant parameters to zero and by setting Me = Mc. For these 'isotropic' analyses, the rotational hardening law (Eq. (2)) is 'switched off and K0 will be overpredicted, analogous to the modified Cam-Clay (MCC) model.

(2) Anisotropy has been 'switched on' while destructuration is 'switched off.

(3) Both the anisotropy and bonding/destructuration parameters have been 'switched on'. The intrinsic compression and creep indices quoted in Table 1 have been used for the analyses incorporating destructuration.

ISO is used hereafter to denote the isotropic response of the Creep-SCLAY1S model with anisotropy and destructuration 'switched off, ANIS is used to denote the anisotropic response with destructuration 'switched off, and A&D denotes the response with both anisotropy and bonding/destructuration 'switched on'.

5. Modelling stone columns using the finite element method

5.1. Previous numerical studies

The majority of numerical studies investigating stone column behaviour have used 2D analysis techniques, e.g. plane strain (Gab et al., 2008) or axisymmetry (Ambily and Gandhi, 2007); 3D modelling has been used by Weber et al. (2008), Kamrat-Pietraszewska and Karstunen (2009), and Killeen and McCabe (2014). In the majority of cases, either the Mohr Coulomb (MC) or HS models have been used to represent the

behaviour of both the granular column material and the treated soil (e.g. Ellouze and Bouassida, 2009; Killeen and McCabe, 2014). Kamrat-Pietraszewska and Karstunen (2009) used the MCC, S-CLAY1, and S-CLAY1S models for the soil and the HS model for the column. Kamrat-Pietraszewska (2011) and Sexton and McCabe (2013) were among the first to use a model incorporating viscous effects with an application to stone columns.

The majority of numerical studies investigating stone column behaviour have declined to use interface elements at the boundary between the granular column material and the in situ soil (e.g. Ambily and Gandhi, 2007; Domingues et al., 2007a, b; Gab et al., 2008), consistent with field observations that columns are tightly interlocked with the surrounding soil due to the lateral displacement caused by re-lowering the poker during column installation (e.g. McCabe et al., 2009). This assumption of full contact at the column-soil interface is adopted in this study.

5.2. Axisymmetric unit cell concept

The axisymmetric unit cell concept is used in this paper, representative of the behaviour of a large grid of regularly-spaced columns subjected to a uniform load, as would be used to support an embankment or a large floor slab, for example. The extent of treatment is usually measured using the area-replacement ratio:

A = I (&" A g I s

where Ac is the cross-sectional area of a single stone column; A is the cross-sectional area of its 'unit cell'; s and Dc denote the column spacing and column diameter, respectively; and g is a constant depending on the column arrangement (Fig. 2). The boundary conditions applied to the unit cell are shown in Fig. 3.

Column (Diameter = Dc, Area

Triangular grid

g = (2V3)/ji

Hexagonal grid g = (3V3)/jt

Fig. 2. Typical column grids encountered in practice. (a) Triangular; (b) Square; (c) Hexagonal.

100 kPa

Upper Carse Clay

Stone Column

0.00 m

-1.50 m -2.50 m

Water table at-1.00 m

® = stress point in column (3)= stress point in soil -8.50 m

L-—____ Fixed in radial direcrion

-14.50m

_ Fixed in vertical and

radial directions

Unit cell radius

Fig. 3. Bothkennar soil profile (not to scale).

5.3. Stone column material parameters

The hyperbolic elastoplastic HS model has been used to model the column material. The parameters have also been derived from Killeen and McCabe (2014) (see Table 3). The model has two yield surfaces: a shear hardening yield surface to incorporate shear hardening and a cap yield surface for compression hardening. The sizes of these yield surfaces are governed by the secant (E50) and oedometric (Eoed) moduli, respectively; the unload-reload modulus, Eur, is used to control the elastic unload-

reload behaviour. The model accounts for the stress dependency of soil stiffness using the following power law:

E = E"

where m is the power dictating the stress dependency of soil stiffness and Eref is a reference stiffness modulus corresponding to a reference pressure,

Table 3. Properties of granular column material.

Y (kN/m3) ф' (°) ¥ (°) c' (kPa)_EOf (MPa)_E5rf (MPa)_Euf (MPa)_p"' (kPa) m_kx (m/d)_k, (m/d)

19 45 15 1 70 70 210 100 0.3 1.7 1.7

6. Axisymmetric 2D modelling

The impact of creep on stone column settlement performance is examined by performing two sets of analyses: one set using the standard soil properties in Table 1 and the other using very low creep coefficients (ц* < 1% of the standard value), in effect eliminating most of the creep effects. It is not possible to use ц = 0 as it would result in division by zero.

These two sets of analyses (referred to subsequently as 'creep' and 'no creep') are carried out with a view to deriving settlement improvement factors (n values) for the three separate scenarios laid out in Section 4.2:

(1) The n values with and without creep for the isotropic case are denoted by nTOTAL(iS0) and npRiMARY(iSo), respectively.

(2) The n values with and without creep for the anisotropic case are denoted by nT0TAL(ANIS) and nPRIMARY(ANIS), respectively. Direct comparison with the isotropic results enables the influence of anisotropy to be established.

(3) The n values with and without creep for the analyses incorporating anisotropy, bonding, and destructuration are denoted by nT0TAL(A&D)

and nPRiMARY(A&D), respectively. The influence of soil destructuration can then be examined.

The variations of radial, vertical, and hoop stress with time and depth

corresponding to scenarios (1)-(3) can be expounded.

The general analysis stages are as follows:

(1) Generate initial stresses using the K0 procedure (Brinkgreve et al., 2011).

(2) Install the stone columns in undrained conditions using the 'wished-in-place' technique (any changes in stresses and state parameters due to column installation are not accounted for). Any out-of-equilibrium stresses generated by the 'wished-in-place' installation are restored using a plastic nil-step.

(3) Apply a load (pa) of 100 kPa in undrained conditions through a plate element placed over the surface of the unit cell. The plate acts as a loading platform and prevents differential settlements at the surface between the column and the soil.

(4) Allow a consolidation phase; settlements effectively cease after full pore pressure dissipation for the 'no creep' case.

7. Computational results and discussion

7.1. Time-settlement behaviour

Settlement versus time plots in logarithmic scale for the untreated 'creep' and 'no creep' cases for the three different scenarios are presented in Fig. 4. The time-settlement behaviour for the isotropic and anisotropic

cases is almost identical. This will be the case for 1D loading if the anisotropy parameters are derived based on the K0 state. The end-of-primary (EOP) consolidation times for the 'no creep' and 'creep' cases are approximately 15,000 days (~40 years) and 40,000 days (~100 years), respectively. The EOP consolidation times are shorter and settlements are lower for the case with bonding and destructuration (since l* << X" and

m << Л

No Creep (ISO)

• No Creep (ANIS)

• No Creep (A&D) Creep (ISO) Creep (ANIS) Creep (A&D)

Time (d) 100 1000

100000

1000000

the isotropic anisotropic ca is almost iden

ent r and ses ical

Fig. 4. Settlement vs. log(time) plots for untreated soil.

The time-settlement behaviour for the treated case (at different reciprocal area-replacement ratios, A/Ac) is compared to that of the untreated case in Fig. 5a and b for the 'no creep' and 'creep' cases, respectively; these plots pertain to the isotropic case. It is evident that the granular columns significantly accelerate the consolidation process; the consolidation time reduces with increasing stone replacement. The findings are consistent with those reported by Kok Shien (2013), who also modelled stone columns using the axisymmetric unit cell concept, albeit using a different soil model and soil profile. Settlement-log(time) plots for the 'anisotropy' and 'anisotropy and destructuration' cases are not presented here because the patterns are relatively consistent with the isotropic case; settlement differences will be reflected in the relevant n values.

7.2. Evolution of settlement improvement factor with time

The evolution of n with time for the 'creep' and 'no creep' cases is plotted in Fig. 6 for the different scenarios; A/Ac = 6 is chosen for illustrative purposes. In all cases, the predicted n values are less than

unity initially because the settlement of treated ground occurs more rapidly than that of untreated ground; however, these n values are of no practical significance. Regardless of the scenario, the 'steady-state' nTOTAL values after EOP are less than the corresponding nPRiMARY values; this holds at all values of A/Ac and is consistent with the findings of Sexton and McCabe (2013, 2015, 2016): when creep is present, settlement improvement factors are lower. 7.3. Comparison of settlement improvement factors

The nPRIMARY and nTOTAL values (after EOP) for the different scenarios are presented in Fig. 7a and b, respectively. At all values of A/Ac, the highest n values arise for the isotropic case and the lowest for the analyses incorporating bonding and destructuration. For the isotropic analyses, K0 is overpredicted (see Section 4.2), resulting in larger horizontal soil stresses which provide more resistance against column bulging. This results in lower settlements for the treated case, and since the settlements of untreated soil for the isotropic and anisotropic cases are similar, npRiMARY(iso) > npRiMARY(ANis).

Time (d) 100 1000

10000 100000 1000000

Time (d)

100000

1000000

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2

Fig. 5. Settlement vs. log(time) plots for the isotropic case. (a) Very low creep coefficient; (b) Standard creep coefficient.

Fig. 6. Evolution of n with time at A/Ac = 6. (a) ISO; (b) ANIS; (c) A&D.

For the analyses incorporating bonding and destructuration (ll << A* and mm << Л the ratio of column stiffness to soil stiffness is lower and slightly lower n values would be expected. However, the n values in Fig. 7a and b are still much lower than those expected because n is not very sensitive to soil stiffness (and hence the ratio of column stiffness to soil stiffness) above a threshold value. The lower n values can be explained as follows:

(1) Increased settlement leads to a reduction in the bonding parameter, X.

(2) Creep (i.e. more settlement) leads to additional bond degradation (for both the untreated and treated cases), and hence a larger reduction in X.

(3) For the 'no creep' case, columns reduce settlement to a larger extent than they do for the 'creep' case (and hence they curtail the amount by which X is reduced).

The extent to which destructuration should be accounted for in design will depend on the initial amount of bonding, X0, which is dictated by the soil sensitivity, St (see Eq. (9)).

Analytical predictions obtained using Priebe (1995), Castro and Sagaseta (2009), and Pulko et al. (2011) are superimposed with the nPRIMARY values in Fig. 7a for comparison. The n values predicted by Castro and Sagaseta (2009) and Pulko et al. (2011) fall between the nPRIMARY(ISo) and nPRIMARY(ANIS) for 4 < A/Ac < 15.

'Creep' settlement improvement factors for the different scenarios are compared in Fig. 7c. These 'creep' settlement improvement factors have been derived based on the slopes of the settlement-log(time) plots after EOP:

nCREEP = ^untreated / ^tre.ted (13)

where mi*ntrelted and ALum denote the slopes of the untreated and treated settlement-log(time) plots, respectively. For each scenario considered, these nCREEP values are lower than the corresponding nPRIMARY values. However, given that the nCREEP values are greater than 1, the columns have a positive impact on reducing long-term creep settlements. The nTOTAL values in Fig. 7b are effectively a weighted average of the nPRIMARY and nCREEP values, dependent on the relative percentages of primary/creep settlement. In general, larger differences between nPRIMARY and nTOTAL would be observed in situations where nPRIMARY is larger to begin with. This occurs because the nPRIMARY values are much greater than the nCREEP values and so a larger effect would be seen in the weighted average.

The relative differences between the nTOTAL and nPRIMARY values for each scenario are investigated in Fig. 8 by plotting (nTOTAL-1) against (nPRIMARY—1) at different values of A/Ac; for an untreated soil, both nTOTAL and nPRIMARY will be equal to 1. Each data point in Fig. 8 corresponds to a nPRIMARY value and a nTOTAL value at a single value of A/Ac. Best-fit lines have been added to each figure, along with their corresponding coefficients of determination (R2). The relationship takes the following form:

ntotal — 1 _ b(nprimary — 1) (14)

where P is the slope of the line. For the isotropic and anisotropic cases, the P values are almost identical, suggesting that the relative values of nTOTAL and nPRIMARY are rather independent of anisotropy. For the case incorporating anisotropy and bonding/destructuration, the value of P is higher because (i) the nPRIMARY(A&D) values are lower to begin with and (ii) m << so the weighted effect of creep is less visible.

/V ///

'•7 i:'

Ч Ü1

i:, i:

tu ¡j

M ►ij

6.0 5.0 4.0

2.0 1.0

4.0 3.5 3.0

s 1.5 1.0

0.5 0.0

3.0 2.5 2.0 1.5 1.0 0.5 0.0

y = R2 0.6225x = 0.9974

"primary

6.0 -1

y = 0.6168x

R2 = 0.9948

y = 0.7691x

R2 = 0.9968

"primary-1

(c) X.J

Fig. 8. ("total 1) vs. ("primary-1 ) plots. (a) ISO; (b) ANIS; (c) A&D.

7.4. Variations of vertical, radial, and hoop stress with time

For the 'no creep' case, the stresses on the soil and column are constant after EOP. For the 'creep' case, vertical stress is transferred from the soil to the column as the soil creeps. This is illustrated in Fig. 9 by plotting the variations of vertical stress ( s'yy ) with time for A/Ac = 3 (for the isotropic case) at points C and S (mid-depth of the lower Carse clay layer) in Fig. 3. The spacing represented by A/Ac = 3 has been selected for presentation purposes because the stress transfer is most pronounced at close spacings, although the same trend holds for the range 3 < A/Ac < 15 considered. The additional stress transferred to the already yielded column results in additional yielding, and hence lower " values for the 'creep' case.

The corresponding variations of radial (s'xx ) and hoop (s'zz) stress in the soil with time are plotted in Fig. 10a and b, respectively. For the 'no creep' case, both s'a and sZ, are constant after EOP. For the 'creep' case, these stresses continue to reduce after EOP, with sZ, reducing to a greater extent than s'a . The reduction in s'xx means that the lateral support imparted onto the column by the soil diminishes due to creep. This leads to additional column bulging, more settlement, and ultimately, a lower load-carrying capacity. The s'a and sZ, reductions will be discussed in more detail in Sections 7.5.2 and 7.5.3.

7.5. Profiles of stress and strain with depth

In this section, distributions of stress and strain with depth in the soil for A/Ac = 3 (with and without creep) are compared after 100 years (after EOp for the untreated case) to highlight the effect of the different features (e.g. anisotropy, bonding and destructuration). The stress and strain profiles have been obtained at the same radius from the column centre as a vertical plane through point S in Fig. 3.

® 7.5.1. Vertical stress profiles

The vertical stress profiles in the soil for A/Ac = 3 are plotted in Fig. 11 for the different scenarios without and with creep, respectively. The stress profiles for the untreated case are also included on the figure to provide a frame of reference for comparison between the 'creep' and 'no creep' cases. The Updated Mesh option (e.g. McMeeking and Rice, 1975) has been used for these analyses and so the final ground surface will be 'lower' when there is more settlement, e.g. for the analyses incorporating creep.

400 350 300 250

iP 200

_ £ b

150 100 50 0

Column

A/Ac=3 (No Creep) at S

A/Ac=3 (Creep) at S

A/Ac=3 (No Creep) at C

A/Ac=3 (Creep) at C

10 100 1000 10000 100000 1000000 Time (d)

Fig. 9. Variation of vertical stress ( s' ) with time at mid-depth for the isotropic case.

70 60 50

ь 30 20 10 0

/ mm — • A/Ac=3 (No Creep) at S

_ — A/Ac=3 (Creep) at S

10 100 1000 10000 100000 1000000 Time (d)

^ — • A/Ac=3 (No Creep) at S

— — A/Ac=3 (Creep) at S

0.1 1 10 100 1000 10000 100000 1000000

Time (d)

Fig. 10. Variations of (a) radial stress ( ) and (b) hoop stress ( o'zz ) in the soil with time at mid-depth.

For each scenario, the columns reduce the vertical stress carried by the soil over the full column length in comparison with the untreated case. For the 'creep' case (Fig. 11b), the stress reductions are larger and increase with depth. These are illustrated using arrows and markers for visual purposes. The surplus stress unloaded from the soil is transferred to the column, as discussed in Section 7.4. The stress transfer (i.e. the stress reduction in the soil) is the smallest for the analyses incorporating destructuration (Fig. 11b) because of the lower creep coefficient (i.e.

m << Ц*).

7.5.2. Radial stress profiles

The corresponding radial stress profiles in the soil for A/Ac = 3 are presented in Fig. 12; reference lines (with no physical significance) have been included on this figure for ease of comparison between the 'creep' and 'no creep' cases. The stress profiles for the untreated case have not been included on these figures for the sake of clarity; the radial stresses in the soil for the untreated case are almost equivalent to those in the soil for the 'no creep' case (Fig. 12a), apart from the isotropic case, for

which the untreated radial stresses are overpredicted. For the 'creep' case (Fig. 12b), the radial stresses in the soil are lower than those for the 'no creep' case (Fig. 12a).

The magnitudes of the radial stress reductions are relatively consistent for the three scenarios considered, although there are significant stress oscillations for the 'creep' case. These oscillations are caused by shear-plane formation in the column (due to additional yielding, e.g. Fig. 13, illustrated for the isotropic case) and extend much deeper for the anisotropic case than for the isotropic case. For the isotropic case, the horizontal stresses in the soil are overpredicted and so additional resistance is provided against lateral column bulging; this inhibits the formation of shear-planes in the column. The magnitudes of the shear-planes are similar for the 'anisotropy and destructuration' case, despite the lower creep coefficient ( m << for this scenario, columns do not have the same beneficial effect on x, as discussed in Section 7.3.

i0 ii i2

i0 ii i2

100 Syy (kPa) (a)

Fig. 11. Profiles of vertical stress in the soil for AtAc = 3. (a) No creep; (b) Creep.

i00 o'y (kPa) (b)

I 7 .g

Si 8 Q

9 i0 ii i2

« i \

........ \

ANIS » \

A&I 3 \

I 7 .g

9 i0 ii i2

\ Ai \

aJ^ 6 1

ISO \ \

ANIS \ \

A&D i'1 \

40 60 o'^ (kPa)

40 60 (kPa) (b)

Fig. 12. Profiles of radial stress in the soil for AtAc = 3. (a) No creep; (b) Creep.

[*MJ] 320.00

i ' Si ? Si ¡> Si1"

Is ■if

I|§ ft?'^ tH™

(a) (b)

Fig. 13. Total shear strains for A/Ac = 3 (isotropic case). (a) No creep; (b) Creep.

The shear-plane formations can also be identified by examining profiles of radial strain (e„) in the soil with depth (Fig. 14). For the 'creep' case, there is a sharp decrease of radial strain at the base of the unit cell (Fig. 14b); this explains the 'jump' in the radial stress profile at the base in Fig. 12b. Additional analyses have shown that these 'jumps' also occur at the base of floating columns; floating columns punch into the underlying soil and so there is also a sudden lateral/radial strain reduction. This is illustrated in Fig. 15 for the isotropic case. 7.5.3. Hoop stress profiles

The corresponding distributions of hoop stress in the soil with depth for A/Ac = 3 are presented in Fig. 16. The plots are presented in a similar format to those for radial stress. The hoop stresses are equal to the radial stresses for the untreated case. The following points are relevant:

(1) The hoop stresses in the soil are lower than the radial stresses (comparing Fig. 16 with Fig. 12) and hence lower than the corresponding hoop stresses in the soil for the untreated case.

(2) The hoop stress reductions for the 'no creep' case are more prominent in the regions where the strains are the largest, e.g. at the depth of the maximum column bulging (which can be approximated from Figs. 13 and 14 as being between 3 m and 4 m below ground level). The hoop stress reductions are caused by plastic deformation which leads to the dissipation of energy; there is more plastic deformation at the bulging depth than that at the base.

(3) The hoop stress reductions are larger for the 'creep' case (Fig. 16b) because there is more plastic deformation throughout the full depth of the profile. A larger hoop stress reduction occurs for the anisotropic case than for the isotropic case (additional plastic deformation and shear plane formation lead to a larger hoop stress reduction, e.g. Fig. 14b). The hoop stress profile for the analyses incorporating anisotropy and bonding/destructuration is approximately parallel to that for the anisotropic case, although the hoop stress reduction is lower (less plastic deformation because m << e.g. Fig. 14b).

..♦•i*

......< ff

. * V ** O*

« . ° o • • A

i «A TPk.....

.......¿1

ISO >

10 ii i2

sxx sxx

(a) (b)

Fig. 14. Profiles of radial strain in the soil for A/Ac = 3. (a) No creep; (b) Creep.

.....v

.....t

........3

Floating

End- bearing

7 1 .g

10 11 12

A A -J

«o* 1 " .

44 A " A- A

4 c>.....

s [.....

** 4 ^ A A A *

Floating

End-bearing """ - -

7 1 .g

10 11 12

Fig. 15. Profiles of radial strain in the soil for A/Ac = 3 (isotropic case): floating vs. end-bearing columns. (a) No creep; (b) Creep.

)mt( h7 tp

9 10 11 12

JL> \ Ya x

ANIS \ i \

A&E \ x \

10 11 12

40 60 c'zz (kPa) (a)

80 100

40 60 c'zz (kPa) (b)

Fig. 16. Profiles of hoop stress in the soil for A/Ac = 3. (a) No creep; (b) Creep.

8. Conclusions and recommendations for future research

A series of axisymmetric analyses has been carried out in conjunction with the elasto-viscoplastic Creep-SCLAY1S model to assess the effectiveness of stone columns in soft creep-prone soils. The columns were wished-in-place and hence installation effects have not been accounted for. Three different scenarios have been considered: (i) isotropy, (ii) anisotropy, and (iii) anisotropy and bonding/destructuration. The main findings are as follows:

(1) For all three scenarios, incorporating creep leads to lower 'total' settlement improvement factors, i.e. lower than 'primary' settlement improvement factors. The 'total' settlement improvement factors are effectively a weighted average of 'primary' and 'creep' settlement improvement factors; the latter are much lower than the former but are, nevertheless, greater than unity. If creep constitutes a significant proportion of total settlement, lower settlement improvement factors for creep settlements should be used in design.

(2) The ratios of 'total' to 'primary' settlement improvement factors are almost identical for the 'isotropic' and 'anisotropic' cases, suggesting the effectiveness of stone columns at arresting creep settlements is independent of anisotropy. A smaller ratio is observed for the case incorporating 'anisotropy and bonding/destructuration' because (i) the 'primary' settlement improvement factors for this scenario are lower to begin with and (ii) the 'intrinsic' creep index is less than the creep index for natural clay (m* << P*) so the effect of creep on the weighted average is less visible.

(3) For the 'creep' case, vertical stress is transferred from the 'creeping' soil to the granular column. The additional vertical stress transferred to the already yielded column causes additional yielding and

explains why 'total' settlement improvement factors are lower than their 'primary' counterparts.

(4) The actual " values (both 'primary' and 'total') are lower for the analyses incorporating bonding and destructuration. The extent to which destructuration should be accounted for in design will depend on the initial sensitivity of the clay; in highly sensitive clays, destructuration will be a greater consideration.

(5) In addition to the vertical stress transfer process, the radial and hoop stresses in the soil also reduce for the 'creep' case. The radial stress reduction results in additional column bulging and a lower load-carrying capacity. The hoop stress reduction is more of a secondary effect, caused by additional plastic deformation (and hence energy dissipation) for the 'creep' case.

(6) The simulations ignored any installation effects, which inevitably are significant (e.g. Castro and Karstunen, 2010; Castro et al., 2014). Hence for future work, it is recommended to complement the present study with analyses that account for installation effects for a more complete understanding of how the stone columns behave in soft sensitive creep-prone soils.

Conflict of interest

We wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.

Acknowledgements

The authors would like to acknowledge the funding provided by the Irish Research Council (IRC) for the research into stone column behaviour in creep-prone soils. The development of the soil model used

herein was carried out as part of CREEP (Creep of Geomaterials, PIAP-GA-2011-286397) project supported by the European Community through the programme Marie Curie Industry-Academia Partnerships and Pathways (IAPP) under the 7th Framework Programme. The support from the BIG (Better Interaction in Geotechnics) project from the Swedish Transport Administration is also gratefully acknowledged.

Notation

The followi"g symbols are used i" this paper:

A = Cross-sectional area of soil unit treated with granular material Ac = Cross-sectional area of granular column Ac/A = Area-replacement ratio Cs = Swelling index Cc = Compression index

Ca = Coefficient of secondary compression / creep coefficient

c£ = Effective cohesion

Dc = Column diameter

E50 = Secant/triaxial modulus

Eoed = Oedometric modulus

Eur = Unload-reload modulus

e0 = Initial void ratio

K0 = Coefficient of lateral earth pressure at rest

K0nc = Coefficient of lateral earth pressure in the normally consolidated condition

g = Constant dependent on column arrangement (square, triangular, or hexagonal)

k, kx, ky = Permeability, horizontal permeability, vertical permeability M, Mc, Me = Slope of CSL, slope of CSL in compression, slope of CSL in extension

M(0) = Stress ratio at critical state

m = Power dictating the stress dependency of soil stiffness (HS model)

" = Settlement improvement factor, " = ¿imMed/^treated

"TOTAL = 'Total' settlement improvement factor (i.e. primary + creep)

"creep = 'Creep' settlement improvement factor

"PRIMARY = 'Primary' settlement improvement factor

"2 = Priebe's (1995) settlement improvement factor

p, p£ = Mean principal total stress, mean principal effective stress

pa = Applied load / load level

pp = Preconsolidation stress / pressure (3D)

pref = Reference pressure

q = Deviatoric stress

Rc = Column radius

St = Sensitivity

s = Column spacing

a0, a = Initial yield surface inclination, yield surface inclination

Y = Bulk unit weight

5 = Settlement

sxx = Radial strain

0 = Lode angle

k, k = Swelling indices

A, A* = Compression indices

At, I* = Intrinsic compression indices

p, p* = Creep coefficients/indices

v = Poisson's ratio

£ = Rate of destructuration

= Effectiveness of shear and volumetric strains in destroying the bonding

s0z = Initial effective stress / pressure (1D) sZ = Preconsolidation stress / pressure (1D)

sXx , s'y, sZ = Effective radial, vertical (axial), and hoop (tangential) stresses

( = Friction angle

rn = Rate of yield surface rotation

md = Effectiveness of shear and volumetric strains in rotating the yield surface

Xo = Initial amount of bonding ^ = Dilatancy angle A = Rate of viscoplastic multiplier

References

Allman MA, Atkinson JH. Mechanical properties of reconstituted Bothkennar soil. Géotechnique 1992; 42(2): 289-301.

Ambily AP, Gandhi SR. Behavior of stone columns based on experimental and FEM analysis. Journal of Geotechnical and Geoenvironmental Engineering 2007; 133(4): 405-15.

Atkinson JH, Allman MA, Boese RJ. Influence of laboratory sample preparation procedures on the strength and stiffness of intact Bothkennar soil recovered using the Laval sampler. Géotechnique 1992; 42(2): 349-54.

Black JA, Sivakumar V, Bell A. The settlement performance of stone column foundations. Géotechnique 2011; 61(11): 909-22.

Brinkgreve RBJ, Swolfs WM, Engin E. PLAXIS 2D 2010. Delft, the Netherlands: PLAXIS B.V., 2011.

Castro J, Karstunen M, Sivasithamparam N, Sagaseta C. Numerical analyses of stone column installation in Bothkennar clay. In: Proceedings of the International Conference on Installation Effects in Geotechnical Engineering (ICIEGE), Rotterdam, the Netherlands, 2013. pp. 212-8.

Castro J, Karstunen M, Sivasithamparam N. Influence of stone column installation on settlement reduction. Computers and Geotechnics 2014; 59: 87-97.

Castro J, Karstunen M. Numerical simulations of stone column installation. Canadian Geotechnical Journal 2010; 47(10): 1127-38.

Castro J, Sagaseta C. Consolidation around stone columns: Influence of column deformation. International Journal of Numerical and Analytical Methods in

Geomechanics 2009; 33(7): 851-77.

Domingues TS, Borges JL, Cardoso AS. Stone columns in embankments on soft soils: Analysis of the effects of the gravel deformability. In: Proceedings of the 14th European Conference on Soil Mechanics and Geotechnical Engineering, Madrid, Spain,

2007a. pp. 1445-50.

Domingues TS, Borges JL, Cardoso AS. Parametric study of stone columns in embankments on soft soils by finite element method. In: Proceedings of the 5th International Workshop on Applications of Computational Mechanics in Geotechnical Engineering, Guimaraes, Portugal, 2007b. pp. 281-91.

Ellouze S, Bouassida M. Prediction of the settlement of reinforced soft clay by a group of stone columns. In: Proceedings of the 2nd International Conference on New Developments in Soil Mechanics and Geotechnical Engineering, Nicosia, North Cyprus, 2009. pp. 182-7.

Gab M, Schweiger HF, Kamrat-Pietraszewska D, Karstunen M. Numerical analysis of a floating stone column foundation using different constitutive models. In: Proceedings of the 2nd International Workshop on the Geotechnics of Soft Soils - Focus on Ground Improvement, Glasgow, 2008. pp. 137-42.

Gens A, Nova R. Conceptual bases for a constitutive model for bonded soils and weak rocks. In: Proceedings of the International Symposium on Geotechnical Engineering of Hard Soils - Soft Rocks, Athens, Greece, 1993, pp. 485-94.

Grimstad G, Degago SA, Nordal S, Karstunen M. Modeling creep and rate effects in structured anisotropic soft clays. Acta Geotechnica 2010; 5(1): 69-81.

Institution of Civil Engineers (ICE). Bothkennar soft clay test site: Characterization and lessons learned. Géotechnique 1992; 42(2): 161-378.

Jardine RJ, Lehane BM, Smith PR, Gildea PA. Vertical loading experiments on rigid pad foundations at Bothkennar. Géotechnique 1995; 45(4): 573-97.

Kamrat-Pietraszewska D, Karstunen M. The behaviour of stone column supported embankment constructed on soft soil. In: Computational Geomechanics, COMGEO I -Proceedings of the 1st International Symposium on Computational Geomechanics, Juan-les-Pins, 2009. pp. 829-41.

Kamrat-Pietraszewska D. Numerical modelling of soft soils improved with stone columns. PhD Thesis. Glasgow, UK: University of Strathclyde, 2011.

Karstunen M, Krenn H, Wheeler SJ, Koskinen M, Zenter R. Effect of anisotropy and déstructuration on the behavior of Murro test embankment. International Journal of Geomechanics 2005; 5(2): 87-97.

Karstunen M, Sivasithamparam N, Brinkgreve RBJ, Bonnier PG. Modelling rate-dependent behaviour of structured clays. In: Proceedings of the International Conference on Installation Effects in Geotechnical Engineering (ICIEGE), Rotterdam, The Netherlands, 2013. pp. 43-50.

Killeen MM, McCabe BA. Settlement performance of pad footings on soft clay supported by stone columns: A numerical study. Soils and Foundations 2014; 54(4): 760-76.

Kok Shien NG. Numerical study of floating stone columns. PhD Thesis. Singapore: National University of Singapore, 2013.

Koskinen M, Karstunen M, Wheeler SJ. Modelling destructuration and anisotropy of a natural soft clay. In: Proceedings of the 5th European Conference on Numerical Methods in Geotechnical Engineering (NUMGE 2002), Paris, 2002. pp. 11 -20.

Leoni M, Karstunen M, Vermeer PA. Anisotropic creep model for soft soils. Géotechnique 2008; 58(3): 215-26.

Leroueil S, Lerat P, Hight DW, Powell JJM. Hydraulic conductivity of a recent estuarine silty clay at Bothkennar. Géotechnique 1992; 42(2): 275-88.

Leroueil S, Tavenas F, Brucy F, La Rochelle P, Roy M. Behaviour of destructured natural clays. Journal of Geotechnical Engineering 1979; 105(6): 759-78.

Liingaard M, Augustesen A, Lade PV. Characterization of models for time-dependent behavior of soils. International Journal of Geomechanics 2004; 4(3): 157-77.

McCabe BA, Nimmons GJ, Egan D. A review of field performance of stone columns in soft soils. Proceedings of the Institution of Civil Engineers - Geotechnical Engineering 2009; 162(6): 323-34.

McMeeking RM, Rice JR. Finite-element formulations for problems of large elastic-plastic deformation. International Journal of Solids and Structures 1975; 11(5): 601-16.

Mesri G. Coefficient of secondary compression. Journal of the Soil Mechanics and Foundations Division 1973; 99(1): 123-37.

Mitchell JK, Kelly R. Addressing some current challenges in ground improvement. Proceedings of the Institution of Civil Engineers - Ground Improvement 2013; 166(3): 127-37.

Moorhead MC. Effectiveness of granular columns for containing settlement of foundations supported on soft clay. PhD Thesis. Belfast, UK: Queen's University, 2013.

Nash DFT, Powell JJM, Lloyd IM. Initial investigations of the soft clay test site at Bothkennar. Géotechnique 1992a; 42(2): 163-81.

Nash DFT, Sills GC, Davison LR. One-dimensional consolidation testing of soft clay from Bothkennar. Géotechnique 1992b; 42(2): 241-56.

Olsson M. On rate-dependency of Gothenburg clay. PhD thesis. Gothenburg, Sweden: Chalmers University of Technology, 2013.

Paul MA, Peacock JD, Wood BF. Engineering geology of the Carse clay at the National Soft Clay Research Site, Bothkennar. Géotechnique 1992; 42(2): 183-98.

Perzyna P. Fundamental problems in viscoplasticity. Advances in Applied Mechanics 1966; 9(2): 244-377.

Perzyna P. The constitutive equations for work-hardening and rate sensitive plastic materials. Proceedings of Vibration Problems 1963; 4: 281-90.

Priebe HJ. The design of vibro replacement. Ground Engineering 1995; 28(10): 31-7.

Pulko B, Majes B, Logar J. Geosynthetic-encased stone columns: Analytical calculation model. Geotextiles and Geomembranes 2011; 29(1): 29-39.

Schanz T, Vermeer PA, Bonnier PG. The hardening soil model: Formulation and verification. In: Beyond 2000 in Computational Geotechnics. Ten Years of PLAXIS. Amsterdam, The Netherlands: A.A. Balkema; 1999. pp. 281-90.

Sexton BG, McCabe BA. Modeling stone column installation in an elasto-viscoplastic soil. International Journal of Geotechnical Engineering 2015; 9(5): 500-12.

Sexton BG, McCabe BA. Numerical modelling of the improvements to primary and creep settlements offered by granular columns. Acta Geotechnica 2013; 8(4): 447-64.

Sexton BG, McCabe BA. Stone column effectiveness in soils with creep: A numerical study. Geomechanics and Geoengineering 2016; Doi:

10.1080/17486025.2016.1151556.

Sexton BG. The influence of creep on the settlement of foundations supported by stone columns. PhD thesis. Galway, Ireland: National University of Ireland, 2014.

Sivasithamparam N, Karstunen M, Bonnier P. Modelling creep behaviour of anisotropic soft soils. Computers and Geotechnics 2015; 69: 46-57.

Sivasithamparam N, Karstunen M, Brinkgreve RBJ, Bonnier PG. Comparison of two anisotropic models at element level. In: Proceedings of the International Conference on Installation Effects in Geotechnical Engineering (ICIEGE), Rotterdam, The Netherlands, 2013. pp. 72-8.

Vermeer PA, Neher HP. A soft soil model that accounts for creep. In: Beyond 2000 in Computational Geotechnics. Ten Years of PLAXIS. Amsterdam, the Netherlands: A.A. Balkema; 1999. pp. 249-61.

Vermeer PA, Stolle DFE, Bonnier PG. From the classical theory of secondary compression to modern creep analysis. In: Proceedings of the 9th International Conference on Computer Methods and Advances in Geomechanics, Wuhan, China, 1998. pp. 246978.

Weber TM, Springman SM, Gab M, Racansky V, Schweiger HF. Numerical modelling of stone columns in soft clay under an embankment. In: Proceedings of the 2nd International Workshop on the Geotechnics of Soft Soils - Focus on Ground Improvement, Glasgow, 2008. pp. 305-11.

Wheeler SJ, Naatanen A, Karstunen M, Lojander M. An anisotropic elastoplastic model for soft clays. Canadian Geotechnical Journal 2003; 40(2): 403-18.

Wheeler SJ. A rotational hardening elasto-plastic model for clays. In: Proceedings of the 14th International Conference on Soil Mechanics and Foundation Engineering, Hamburg, 1997. pp. 431-4.

Yin JH, Zhu JG, Graham J. A new elastic viscoplastic model for time-dependent behaviour of normally and overconsolidated clays: Theory and verification. Canadian Geotechnical Journal 2002; 39(1): 157-73.

Yin ZY, Chang CS, Karstunen M, Hicher PY. An anisotropic elastic-viscoplastic model for soft clays. International Journal of Solids and Structures 2010; 47(5): 665-77.

Yin ZY, Karstunen M, Chang CS, Koskinen M, Lojander M. Modeling time-dependent behavior of soft sensitive clay. Journal of Geotechnical and Geoenvironmental Engineering 2011; 137(11): 1103-13.

Yin ZY, Karstunen M. Influence of anisotropy, destructuration and viscosity on the behavior of an embankment on soft clay. In: Proceedings of the 12th International Conference of the International Association for Computer Methods and Advances in Geomechanics (IACMAG), Goa, India, 2008. pp. 4728-35.

Zentar R, Karstunen M, Wiltafafsky C, Schweiger HF, Koskinen M. Comparison of two approaches for modelling anisotropy of soft clays. In: Proceedings of the 8th International Symposium of Numerical Models in Geomechanics (NUMOG VIII), Rome, 2002. pp. 115-21.

Zhou C, Yin JH, Zhu JG, Cheng CM. Elastic anisotropic viscoplastic modeling of the strain-rate-dependent stress-strain behavior of ^-consolidated natural marine clays in triaxial shear tests. International Journal of Geomechanics 2005; 5(3): 218-32.

Brian Sexton graduated with a BE in Civil Engineering from National University of Ireland (NUI) Galway in October 2010 and then a PhD in Civil Engineering, also at NUI Galway, in October 2014. He is currently working as a Geotechnical Engineer with AGL Consulting in Dublin, Ireland on a wide range of design projects. His research interests include ground improvement techniques, numerical modelling, soft soil behaviour, and retaining wall design. Brian represented Ireland at the 25th European Young Geotechnical Engineers Conference in Sibiu, Romania, in June 2016.

Bryan McCabe completed both his Bachelor's Degree in Civil Structural and Environmental Engineering and PhD in Geotechnical Engineering at Trinity College Dublin, Ireland. Bryan has been a lecturer in Civil Engineering since 2001, with the exception of a one-year secondment to Keller Ground Engineering, Coventry, UK. Bryan is currently Head of Civil Engineering at NUI Galway. His research interests include piling, ground improvement, tunnelling and the Irish pyrite problem.

Minna Karstunen is Professor in Geotechnical Engineering at Chalmers University of Technology, a Fellow of the Institution of Civil Engineers and a Fellow of the Royal Society of Sciences and Arts in Gothenburg. After her PhD research at the University of Wales Swansea, she had academic appointments in e.g. University of Glasgow and Strathclyde. Minna's research combines Applied Mechanics with Geotechnics, motivated by the society's need to design, construct and maintain infrastructure and buildings on/with materials with geological origin, in the most safe, economic and sustainable way. Minna's field of expertise is constitutive and numerical modelling of geomaterials, and she has coordinated numerous EC funded projects related to her research field.

Nallathamby Sivasithamparam is a project engineer and researcher in the field of Computational Geomechanics at Norwegian Geotechnical Institute (NGI), Norway. He received his BSc Eng. (Hons) in Civil Engineering from University of Peradeniya, Sri Lanka in 2004, MEng in Earthquake Engineering from University of Tokyo, Japan in 2007, and PhD in Computational Geomechanics from University of Strathclyde, Scotland in 2012. He was a Marie-Curie Research Fellow at PLAXIS (The Netherlands), National Autonomous University of Mexico (Mexico) and Chalmers University of Technology (Sweden). His main research interests cover development of constitutive models for natural soils and their numerical implementation into finite element codes.