An analytical solution for the settlement of stone columns beneath rigid footings

This paper presents a new approximate solution to study the settlement of rigid footings resting on a soft soil improved with groups of stone columns. The solution development is fully analytical, but finite element analyses are used to verify the validity of some assumptions, such as a simplified geometric model, load distribution with depth and boundary conditions. Groups of stone columns are converted to equivalent single columns with the same cross-sectional area. So, the problem becomes axially symmetric. Soft soil is assumed as linear elastic, but plastic strains are considered in the column using the Mohr–Coulomb yield criterion and a non-associated flow rule, with a constant dilatancy angle. Soil profile is divided into independent horizontal slices, and equilibrium of stresses and compatibility of deformations are imposed in the vertical and horizontal directions. The solution is presented in a closed form and may be easily implemented in a spreadsheet. Comparisons of the proposed solution with numerical analyses show a good agreement for the whole range of common values, which confirms the validity of the solution and its hypotheses. The solution also compares well with a small-scale laboratory test available in the literature.

Abstract This paper presents a new approximate solution to study the settlement of rigid footings resting on a soft soil improved with groups of stone columns. The solution development is fully analytical, but finite element analyses are used to verify the validity of some assumptions, such as a simplified geometric model, load distribution with depth and boundary conditions. Groups of stone columns are converted to equivalent single columns with the same cross-sectional area. So, the problem becomes axially symmetric. Soft soil is assumed as linear elastic, but plastic strains are considered in the column using the Mohr-Coulomb yield criterion and a non-associated flow rule, with a constant dilatancy angle. Soil profile is divided into independent horizontal slices, and equilibrium of stresses and compatibility of deformations are imposed in the vertical and horizontal directions. The solution is presented in a closed form and may be easily implemented in a spreadsheet. Comparisons of the proposed solution with numerical analyses show a good agreement for the whole range of common values, which confirms the validity of the solution and its hypotheses.
The solution also compares well with a small-scale laboratory test available in the literature. Compressive stresses and strains are assumed as positive throughout the paper. Effective stresses are used throughout the paper, which are equal to total stresses but for hydrostatic pore pressures because drained conditions are assumed. For the sake of brevity, dash notation is not used. Stone columns, also known as granular piles or aggregate piers, are commonly employed to improve soft soils for foundation of embankments or structures. They are vertical boreholes in the ground, filled upwards with gravel compacted by means of a vibrator. The inclusion of gravel, which has a higher strength, stiffness and permeability than the natural soft soil, improves the bearing capacity and the stability of embankments and natural slopes, reduces total and differential settlements, accelerates soil consolidation and reduces the liquefaction potential.
Stone columns are usually installed in large uniform grids under embankments or uniform distributed loads. In those cases, taking advantage of symmetries, theoretical studies usually simplified the problem to a ''unit cell'', i.e. only one column and the corresponding surrounding soil. So, there are analytical solutions that provide the settlement of a unit cell [3,18] and its evolution with time [10]. On the other hand, when focus is, for example, on stability of lateral slopes of an embankment, plane strain models, converting the columns to granular trenches, are usually considered [4].
More recently, stone columns have also been deployed beneath small isolated pad or strip foundations at low or moderate loading conditions [25]. The bearing capacity of those groups of stone columns has predominantly been the focus of previous studies [4,14,26,27]. However, stone columns are installed in soft soils that can undergo large displacements at relatively low loads and the serviceability limit state may be critical for their design. Black et al. [5] identified the lack of information regarding the settlement performance and conducted small-scale laboratory tests. Similarly, Shahu and Reddy [21] performed laboratory tests and their numerical simulation. Killeen [12] and Ng [16] have studied groups of stone columns using two-and three-dimensional (2D and 3D) finite element analyses. They studied the influence of several geometric and material properties, such as column length, friction angle, diameter, spacing and stiffness. Recently, the author [8] has also studied numerically groups of stone columns under rigid footings to investigate the influence of the number of columns, column arrangement and column length.
Analytical solutions are useful tools for design, but they are very limited for groups of stone columns beneath rigid footings. Priebe [17] suggests calculating the settlement of a group of columns applying a correction factor, lower than 1, to the settlement of an infinite grid of stone columns, i.e. the settlement of the ''unit cell''. This correction factor is tabulated and depends on the number of columns and the depthto-diameter ratio. Calculation of the tabulated correction factor is based on lower lateral confinement of the outer columns of the column group below the rigid footing, but no further details are provided. Stuedlein and Holtz [23] summarize some simplified solutions, such as those based on spring analogy using the modulus of subgrade reaction [20], and develop a multilinear regression equation, using a database of 58 field load tests to fit the linear coefficients. Ng [16] also proposes a simplified equation fitting coefficients. In this later case, the fitting is done using numerical results.
The review of the analytical solutions to predict the settlement of groups of stone columns beneath rigid footings presented above shows the current interest in the topic and the semi-empirical nature of the proposed solutions, e.g. based on the fitting of experimental [23] or numerical results [16]. Therefore, this paper tries to contribute to more rigorous solutions and presents a new solution that is analytically developed. The analytical treatment of a group of several columns is very complex and some simplifying assumptions are needed, e.g. a simplified geometric model and a simplified stress distribution beneath the footing are used. Those simplifying assumptions are presented in Sect. 2. The new analytical solution is developed in detail in Sect. 3, and a calculation example is presented in Sect. 4. The paper concludes with a validation against numerical results and published experimental data (Sect. 5) and the conclusions.

Model assumptions
The problem to be analysed in this paper is the settlement of a group of stone columns beneath a rigid footing (Fig. 1a). Drained conditions are assumed, so the soil and column responses are studied in effective stresses. However, dashes are not added to stress notation for the sake of brevity. Perfect bonding is assumed between soil and column at their interface, as it is a common practice (e.g. [1]) because stone columns are tightly interlocked with the surrounding soil. Column installation [9] is not considered. A key parameter of the problem is the area replacement ratio, a r , which is defined as the area of the columns, A c , divided by the loaded area, A l . Sometimes, the area replacement ratio is specifically called the footprint replacement ratio to clarify that, at the surface, the loaded area is the area of the footing.

Column length
Columns are assumed to reach a rigid substratum (endbearing columns) or, for floating columns, to have a length, L, equal or higher than the critical column length. Ng [16] and Castro [8] have shown that for a non-stratified soil profile, there is a critical column length, whose value is around 1.5-2 times the footing width, B, for settlement reduction. Castro [8] also shows that for failure conditions, i.e. for bearing pressure, the critical length is lower, around 0.5B. The existence of critical column lengths and their values are justified by the depth of the pressure bulb for settlement reduction and the depth of the failure mechanism for bearing pressure.
The presented analytical model does not consider column punching or penetration into the underlying soil, but that is not possible for end-bearing columns and is small for columns lengths equal or higher than the critical length [8], which will be assumed on the safe side as 2B in this paper.

Simplified geometric model
The analysis of a group of several stone columns in the actual geometry is very complex; therefore, a simplified geometric model is used here. Ng [16] and Castro [8] show that the number of columns and their configuration in a group of stone columns beneath a rigid footing, if a r is kept constant, have little influence on the load-settlement response. This is visible from numerical [8,16] and laboratory results [5,27] and is an important feature because it justifies simplified geometric models that have long been used, such as converting stone columns in concentric rings with the same cross-sectional area, e.g. [15], and homogenization techniques, i.e. an equivalent soil with improved properties, e.g. [19]. Those simplified geometries allow for a 2D axially symmetric analysis (Fig. 1). To get a 2D model, it may also be necessary to transform the footing, e.g. a square footing, to a circular footing with the same area. The transformation of the footing follows a similar principle of same area and is also common practice.
Another simplified 2D geometric model that also neglects the small influence of the number of columns and their configuration is converting all the columns to only one central column with the same cross-sectional area [8] (Fig. 1d). This is the model used here to develop the analytical solution.

Depth and stress distribution
Similar to hand calculations of footings, the soil profile with depth is divided into horizontal slices and a simplified stress distribution is considered. This means that an average constant stress, p a (z), is applied on each horizontal slice at depth, z (Fig. 2). The applied stress on the footing, p a (0), spreads out at a slope, whose value is usually assumed for an unimproved soil as 2V (vertical) to 1H (horizontal). However, stone columns are stiffer than the natural soft soil and concentrate the load, so the load spreads out on a narrower area. Shen and Blackburn [20] propose a slope of 4V:1H from the base of the footing (Fig. 2). So, the average vertical stress at each depth may be obtained as where D is the diameter of the footing. This stress distribution is directly taken from that used for pile groups (e.g. [24]). It is worth clarifying that the spreading of the vertical stress with depth in an improved soil obviously depends on the relative stiffness between the columns and the soil. So, for stone columns, the load is expected to spread on a slightly wider area than that for piles, but the approximate proposal by Shen and Blackburn [20] will be here used as a first approach on the safe side. Additionally, a constant slope is assumed because a homogeneous soil is considered, but for stratified soils, the load distribution slope may vary. For example, if there is an upper crust, usually much stiffer than the soft soil below, the load distribution slope should be higher in that upper crust layer.
On the other hand, the load on the footing displaces the soil downwards, and as the soil is not laterally confined, there is also an important horizontal displacement radially outwards. So, the columns and the soil beneath the footing expand radially (e r \ 0), while the soil far from the footing presents compressive radial strains (e r [ 0). The limit between compressive (positive) and extensive (negative) radial or horizontal strains is shown in Fig. 3 for several cases, which were numerically simulated in Castro [8]. Figure 3 shows that the boundary between compressive and extensive horizontal strains is approximately at a slope of 4V:1H from the base of the footing, which happens to be the same as the assumed vertical stress distribution slope.
The footing is also assumed to be at the ground level, yet it will be usually embedded into the ground, which may be considered as an extra safety or later modified by an embedment correction factor.

Boundary conditions and material properties
The solution is developed for a horizontal slice at a depth z (Fig. 4), and shear stresses between slices at different depths are neglected. The overall behaviour is obtained by means of addition of the solution, i.e. vertical deformation, at different depths. An average vertical stress, p a (z), is applied on each slice and is distributed between soil and column. Average values of the vertical stress on the soil, r zs , and on the column, r zc , are used in the development of the analytical solution. The vertical strain at each horizontal slice, e z , is assumed as uniform, which is an acceptable approximation because the load is applied on a rigid footing.
The boundary condition on the left-hand side ( Fig. 4) is the axis of symmetry, and on the right-hand side, the radial strain is assumed as zero (e r = 0), as explained above and shown in Fig. 3. It is worth noting that this boundary condition (e r = 0) allows for a horizontal displacement radially outwards, which means that there is not a full lateral confinement and the vertical strain is higher than for the ''unit cell'' case, where the lateral displacement at the outer boundary is restrained (e.g. Pulko et al. [18]). The radial distance, at which this boundary condition has to be applied, coincides with the lateral extension of the loaded area, r l (z), which varies with depth as explained in Fig. 2. This means that the area replacement ratio decreases with depth, a r (z), because the area of the column, A c , is constant and the loaded area, A l , increases with depth. The area of the column is A c ¼ Np=4d 2 c ¼ p=4d 2 c;eq , where N is the number of columns in the group, and the subscript ''eq'' refers to the central equivalent column. So, similar to Eq. (1), the area replacement ratio is given by Soil is assumed as elastic (E s , m s ), but plastic strains are also considered in the column using the Mohr-Coulomb yield criterion and a non-associated flow rule, with a constant dilatancy angle (E c , m c , / c , w c ).
The accuracy of the above simplifying assumptions will be evaluated in Sect. 5, comparing the results of the analytical model with numerical analyses. slice at a depth z of the unit cell is analysed (Fig. 4). A horizontal slice of the column is a vertical solid cylinder subjected to a vertical uniform pressure r zc and a radial pressure r rc at its lateral wall; i.e. it is subjected to a triaxial stress state. Its elastic solution is (e.g. [3]) Radial and hoop (circumferential) stresses are the same and constant within the column. Therefore, e rc = e hc = -s r /r c , where s r is the radial displacement at the soil/column interface, which is positive radially outwards. Compressive stresses and strains are considered positive, and Lamé's elastic constants is the shear modulus or second Lamé's constant).
The surrounding soil is a cylinder with a central cylindrical cavity, subjected to a vertical uniform pressure r zs , a radial pressure r rs at the cavity wall (soil/column interface) and a null radial strain (e r = 0) at the external boundary (r = r l ). Its solution is derived in detail in Appendix 1 and is equal to The stresses in Eqs. (3) and (4) correspond to stress increments produced by the applied load, p a . Compatibility of deformation at the soil/column interface in the horizontal direction gives e rc = e hc = -e hs . Imposing radial equilibrium at the soil/column interface (r rc ¼ r rs ), the relationship between horizontal and vertical strains is obtained where To get the vertical strain, the condition of vertical equilibrium to soil and column is imposed.
Then, the vertical strain may be expressed as a function of an equivalent elastic modulus where the equivalent elastic modulus is and E m is the confined (oedometric) elastic modulus. The first and second terms in Eq. (9) are the weighted values of the soil and column moduli, respectively. Vertical and horizontal stresses can be determined substituting Eqs. (8) and (5) into Eq. (3) for the column and into Eq. (4) for the soil.

Plastic column
A key issue to correctly capture stone column behaviour is considering its yielding, which can be adequately modelled with the Mohr-Coulomb yielding criterion and a nonassociated flow rule for the plastic strains, with a constant dilatancy angle w c 6 ¼ u c ð Þ : Stone column material is usually dense and granular with a negligible cohesion.
As for the elastic case, the column is assumed to be in compressive triaxial condition, and then, major principal stress direction is vertical. Now, column strains have an elastic part (Eq. 3) and another plastic (Eq. 11). Superscript ''p'' is used here to denote the plastic component. Combining those two components, the stress-strain relationship is: where C E is the following constant of the column, which has units of stiffness: The column constant C E gives an idea of the importance of the elastic strains in the column during its plastic deformation [11]. If the elastic strains are disregarded, 1/ C E = 0. Now, in the elastic-plastic analysis, the effective stresses in the yield condition (Eq. 10) have to include also the initial stresses, existing before load application. Assuming an initial geostatic stress state, given by the effective unit weight of soil and column, c 0 s and c 0 c , and the coefficient of lateral earth pressure at rest of the soil, K 0s , the initial effective stresses are Subscript ''i'' is used here to denote initial (prior to loading) stresses. Initially, the column is in an elastic state, because K 0s [ K ac always, but the applied load, p a , may cause column yielding. The value of the applied load that causes column yielding, p y a , may be obtained using Eq. (10) r rc;y r zc;y ¼ K 0s c 0 s z þ Dr rc;y c 0 c z þ Dr zc;y where the subscript ''y'' refers to column yielding, and the stress increments (denoted as Dr) are due to load application and may be obtained using Eqs. (3), (5) and (8). Substituting those equations in Eq. (15), p y a is obtained. where So, the applied load may be divided into the yielding load and the plastic increment load, p a ¼ p y a þ p p a . The solution for the plastic increment, p p a , may be developed following the same procedure as for the elastic case, but replacing Eq. (3) by Eq. (12) to account for column elasto-plastic behaviour. Now, the relationship between horizontal and vertical strains is where Imposing the condition of vertical equilibrium to soil and column, the vertical strain is where E p ml is an equivalent modulus for the plastic increment where As for the elastic case, the vertical stresses and the radial displacement at the soil/column interface may be obtained substituting Eqs. (20) and (18) in Eqs. (4) and (12).

Solution
The final solution for the vertical strain, e z , of a horizontal slice at depth z may be expressed as To get the footing settlement, s z , the soil profile beneath the footing has to be divided into several slices of finite thickness, Dt. The vertical deformation of each slice is approximated by the vertical strain at the centre of each slice. The addition of the vertical deformation of all of the slices gives the surface settlement To clarify the calculation procedure, an example is presented in the next section. The solution may be easily implemented in a spreadsheet or any programming language; an example is shown in Appendix 2.

Calculation example
In this example, the soil profile is simplified to only one homogeneous soil layer of thickness H = 10 m. The elastic properties of soil and column are Poisson's ratios of m s = m c = 0.33 and Young's moduli of E s = 2 MPa and E c = 30 MPa, respectively. This means a modular ratio of 15. The friction and dilatancy angles for the column material are assumed as / c = 45°and w c = 10°, respectively. For the sake of simplicity, the ground water level is at the surface and the unit weights of soil and column are c s = c c = 20 kN/m 3 , which means effective unit weights of c 0 s = c 0 c = 10 kN/m 3 . The coefficient of lateral earth pressure at rest of the soil is set equal to K 0s = 0.6, disregarding installation effects [9].
The foundation to be studied is a rough rigid square footing of side length B = 5 m under a vertical pressure of p a = 50 kPa. The soil beneath the footing is improved with 4 stone columns, whose diameter is d c = 0.9 m. So, the total column area is A c = 2 m 2 , and divided by the footing area, A l = 25 m 2 , gives the area replacement ratio at the base of the footing a r (0) = 10 %. As shown in [8,16], the number of columns is not relevant if a r is kept constant.
The square footing (B = 5 m) is converted to a circular footing of same area. Therefore, its diameter is In some practical situations, the diameter of the equivalent circular footing is directly taken as the side length of the square footing. The four columns are also converted to a central equivalent column. Its diameter gives the same total column area, If the columns, and consequently the equivalent central column, are assumed to reach a rigid substratum (end-bearing columns), their length is L = H = 10 m.
As an example, the soil and the equivalent column are divided into five horizontal slices of same thickness, Dt = 2 m. For this case, dividing the soil into five slices is accurate enough. No further comments are provided on the necessary number of horizontal slices as it is similar to other settlement calculations and depends on the natural layering of the soil. The applied vertical stress and the area replacement ratio at the centre of each horizontal slice are calculated using Eqs. (1) and (2). The results are shown in Table 1. Equations (9) and (21) are used to get the equivalent moduli for the elastic and plastic cases, E e ml and E p ml , respectively. They are slightly different for each horizontal slice because they depend on the area replacement ratio. Next, the yielding load, p y a , is obtained with Eq. (16). The vertical strain at the centre of each horizontal slice, e z , is given by Eq. (23). Finally, multiplying the slice thickness (2 m) by the vertical strains gives the deformation of each slice, s zi , and the addition of the deformation of the 5 slices is the footing settlement, s z = 70 mm (Table 1).

Validation and parametric analyses
Using the previous calculation example as a reference case, parametric studies were carried out varying several properties. The calculation example and their parametric variations were also compared with numerical analyses to validate the accuracy of the analytical solution. The finite element method was used for the numerical simulations. Similar to Castro [8], the codes Plaxis 3D 2012 [6] and Plaxis 2D 2012 [7] were used for the full 3D models and the simplified 2D axisymmetric models, respectively. The numerical models for the calculation example are shown in Fig. 5. The results of the full 3D and the simplified 2D models are very similar [8]. So, for the sake of simplicity, only the results of the simplified 2D model are shown unless otherwise stated. The soil and column properties in the numerical model are the same as those used for the analytical solution.

Calculation example
The footing settlement calculated in the previous section (s z = 70 mm) is compared with the footing settlement computed numerically, 59 and 57 mm, using the 2D and 3D models, respectively. The small difference between the 2D and 3D models is at least partly caused by the higher order of the elements used in the 2D model [8]. Although the analytical solution slightly overpredicts the settlement as explained below, it captures well the variation with depth (Fig. 6). In the numerical models, the vertical line cross section corresponds to the centre of the footing, which means that it is in the column for the 2D model and in the soil for the 3D model. The oscillation at 2 m depth for the 2D model is caused by a shear band in the column (e.g. Pulko et al. [18]). To further analyse the calculation example, the horizontal displacements at the soil/column interface are compared in Fig. 7. The comparison is only possible with the 2D model because in the 3D model, there are four columns instead of the central equivalent column and the results are not comparable. Above roughly z = 5 m, the horizontal displacements are significant, while below that depth, they are much lower. This is caused by the plastic strains in the column in the upper 5 m, while in the lower part, the column remains elastic. The analytical solution properly reproduces the horizontal displacements but for the upper (z = 0) and lower (z = 10 m) boundaries, which are rough (s rc = 0) in the numerical model. Those differences in the boundary conditions explain part of the difference in the settlement. Numerical simulations in the 2D model were performed using smooth upper and lower boundaries, and the settlement increased from s z = 59 to 61 mm. So, there is still a difference between the numerical model and the analytical solution, which is mainly caused by the stress distribution assumed for the analytical solution (4V:1H) (Fig. 2), as shown in the following.
The vertical stresses, including the initial geostatic stress, are plotted at several depths in Fig. 8. The assumed Depth, z (m) Fig. 6 Vertical displacement with depth distribution of the vertical stress (Fig. 2) is obviously a simplification, and in the numerical analyses, there is not an abrupt distinction between loaded and not-loaded soil. However, the stress distribution between soil and column is well captured by the analytical solution, and the general agreement is good. For this calculation example, the stress distribution assumed in the analytical solution (4V:1H) slightly overpredicts the vertical stresses, and consequently the vertical strains, particularly for z [ B (e.g. Fig. 8c). A 3V:1H distribution of the vertical load gives a better approximation and matches nearly perfectly the surface settlement (Fig. 9). It is worth noting that a 4H:1V distribution is kept for the external boundary condition (e r = 0). In the following, as an estimation on the safe side, the 4H:1V distribution for both the vertical load and the external boundary condition is used if not otherwise stated.

Influence of soil yielding
The analytical solution assumes an elastic behaviour for the soil. For spread loadings and the unit cell model, the soil reaches its maximum strength only near the column and for high load levels. Consequently, assuming an elastic behaviour for the soil gives good settlement predictions for the unit cell model (e.g. [18]). However, the influence of soil plastic deformations for small groups of columns is more significant (e.g. [8]). To evaluate that influence, numerical analyses were performed. As an example, Fig. 10 shows the load-settlement curves of the calculation example but including the numerical simulation with an elasto-plastic behaviour of the soil. The Mohr-Coulomb yield criterion and a non-associated flow rule, with a constant dilatancy angle, were used for both soil and column. A low strength (/ s = 23°, c s = 3 kPa) and a nondilatant behaviour (w s = 0) were chosen for the soil.
Obviously, plastic strains in the soil increase the settlement but for low vertical loads (p a \ 20 kPa for the calculation example). For high loads, plastic strains in the soil notably increase the settlement; ultimately, the bearing pressure is reached. Nevertheless, the proposed analytical solution is just for serviceability limit states, where soil yielding should be limited. For the calculation example and a service pressure around p a = 50 kPa, the settlement overestimation provided by the assumption of a 4V:1H distribution of the vertical load is compensated by the increase in the settlement that is caused by the soil plastic strains (Fig. 10).

Floating columns
So far, only end-bearing columns have been considered, but the analytical solution is also applicable to floating columns, whose length is equal or higher than critical (L C 2B). Using the calculation example as a starting point, the soft soil layer thickness, H, was increased from 10 to 40 m, while keeping the column length equal to the critical length (L = 10 m). The case with columns (improved soil) and the case without columns (unimproved soil) were simulated, and, in both cases, the footing settlement (s z and s z0 , respectively) increases with H (Fig. 11). The settlement increase beyond H = 40 m = 8B is small. For example, the settlement for a homogeneous elastic half-space without columns is around 96 mm. The settlement reduction factor achieved with the improvement, which is defined as the coefficient between the settlement with columns and without columns, b = s z /s z0 , does not notably change beyond the critical length (Fig. 11) [8,16]. Therefore, the settlement reduction factor of floating columns whose length is equal or higher than critical (L C 2B) may be estimated studying those columns as if they were end-bearing columns. To get the settlement reduction factor, the settlement with columns is obtained using the proposed solution and the settlement without columns may be estimated also analytically using any conventional  Fig. 11 Influence of soil layer thickness. Numerical analyses method (e.g. [13]). In some cases, assuming L C 2B may not be realistic, and then, the proposed solution is not applicable.

Influence of area replacement ratio
In the calculation example, the area replacement ratio at the base of the footing is low, namely a r (0) = 10 %, and therefore, the footing settlement is reduced by a small amount (Fig. 11). If a further reduction in the settlement is required, a r may be increased. Figure 12 compares the footing settlement for several values of a r . The settlement predicted by the proposed analytical solution assuming a 3V:1H load distribution matches very well the settlement computed numerically, while the settlement predicted by the proposed analytical solution assuming a 4V:1H load distribution gives a good estimation on the safe side.

Influence of footing width
A parametric analysis was also performed varying the footing width (Fig. 13). For an infinitely small footing, the settlement is null and it increases with the footing width as the loaded area increases. However, as the factor D/ H increases, the lateral confinement of the columns also increases and the settlement tends to an asymptotic value that corresponds to the ''unit cell'' model (full lateral confinement). As for previous cases, the analytical solution gives a conservative estimation of the computed settlement for a 4V:1H load distribution (Fig. 13a), while a good matching is achieved using a 3V:1H distribution (Fig. 13b).
For high values of D/H, e.g. higher than 1 for a r = 30 %, the lateral confinement is high and the agreement between the proposed solution and the numerical analyses is worse because the hypothesis about the location of the external boundary where e r = 0 is no longer valid. Under high lateral confinement, i.e. high values of D/H, there is not a spreading of the external boundary and that boundary is nearly a vertical line from the edge of the footing (e.g. Fig. 14). This means that the proposed solution has been developed for columns beneath footings or pads, and consequently, the solution and its hypotheses are valid for those cases but not for rafts or large loaded areas.

Influence of soil and column properties
To complete the parametric study and verify that the proposed analytical solution is able to accurately predict the footing settlement for the common range of problem parameters, the soil and column properties were also varied. Firstly, the column strength, namely the friction angle, was varied from 35°up to 50°. The dilatancy angle of the column was also varied accordingly, from 0°up to 15°. The footing settlement evidently decreases as the column friction and dilatancy angles increase (Fig. 15), because the column yields later and the plastic strains are less important. The analytical solution assuming a 4V:1H load distribution gives a good estimation on the safe side; only for high area replacement ratios and very low column strengths (e.g. a r = 50 %, / c = 35°, w c = 0), the analytical solution gives slightly higher settlement. Secondly, the soil and column Young's moduli were varied. For the same modular ratio, E c /E s , the results of both parametric analyses are the same if the footing settlement is expressed in terms of the settlement reduction factor, b (Fig. 16). The results confirm that the proposed analytical solution gives a good estimation on the safe side. The agreement between the analytical solution and the numerical analyses is better for high values of a r and E c /E s , because the load distribution is closer to 4V:1H. On the contrary, for low a r or E c /E s , the vertical load spreads out on a wider area and the load distribution slope of 4V:1H assumed for the analytical solution gives conservative values, i.e. overestimates b. Ultimately, for a r = 0 % or E c /E s = 1, the analytical solution would significantly overestimate b, if a 4H:1V load distribution slope is assumed instead of a more realistic slope of 2H:1V for the case of no improvement.

Comparison with published experimental measurements
To fully validate the proposed analytical solution, a comparison with experimental measurements is here presented. The set of high-quality laboratory experiments performed by Sivakumar et al. [22] were used because of the detailed information provided. Sivakumar Fig. 16 Influence of soil and column stiffnesses column and subjected to foundation loading under drained conditions. The resulting area replacement ratios are a r = 44, 69 and 100 %, respectively. Speswhite kaolin and crushed basalt were used to model the soft soil and the stone column. Their properties may be found in Black et al. [5]. Those relevant for the analytical solution are E s = 4 MPa, E c = 30 MPa, / c = 43°. Common values were assumed for other parameters that are not provided, m c = m s = 0.33 and w c = 10°. The soil samples were isotropically consolidated (K 0 = 1) to 50 kPa. This is the initial stress state for this case (Eq. 14). The soil and column weight were neglected. The columns were fully penetrating (L = 400 mm), but the deformation below 2D = 120 mm was assumed to be small, and therefore, for the analytical solution, the column length was input as 120 mm.
Although there is some scatter in the laboratory data, the predicted settlement by the analytical solution matches reasonably well the laboratory measurements for design footing pressures (Fig. 17). For higher pressures, e.g. p a [ 60-100 kPa, the experimental load-displacement curves bend towards the bearing pressure because of soil plastic deformation.

Limit of the applied load
As shown in Sect. 5.2 and Fig. 17, the proposed solution is only valid for low or moderate applied loads that do not generate significant plastic strains in the natural soft soil. An exhaustive analysis of the soil plastic strains is beyond the scope of the paper, but further information and some comments are presented here to provide some guidance on the limit of the applied load when using the proposed solution.
For the calculation example, the proposed solution gives conservative settlements for loads not exceeding 50 kPa (Fig. 10). This limit increases with the soil strength, footing diameter and area replacement ratio. As an example, Fig. 18 shows the settlement predicted by the analytical solution and the numerical analysis, considering soil plastic strains (/ s = 23°, c s = 3 kPa) for different area replacement ratios. When the predicted values by the analytical solution and the numerical analysis coincide, this is considered the limit of the applied load to be used in the analytical solution. The limit load for the calculation example, i.e. a r = 10 %, is 50 kPa, but this value increases up to 75 kPa for a r = 60 % and up to 100 kPa for a r = 80 %. Those values are higher for higher soil strengths and footing diameters.
In comparison with the laboratory measurements presented in the previous section (Fig. 17), the limit load is between 60 and 100 kPa, depending on a r and the scatter of the laboratory results. Although the footing diameter is small, the limit load is similar to the previous calculation example because of the higher soil strength.
In summary, for common situations (footing diameter, soil strength and a r ), the limit of the applied load to be used in the analytical solution is between 50 and 100 kPa.

Conclusions
A new approximate solution has been developed to study the settlement of rigid footings resting on a soft soil improved by groups of stone columns. The solution development is fully analytical, imposing equilibrium of stresses and compatibility of deformations. Some simplifying assumptions are necessary, such as: • Groups of columns are converted to equivalent single columns with the same cross-sectional area. So, the problem becomes axially symmetric. • Applied load distributes with depth at a 4V:1H slope.
For some cases, a 3V:1H load distribution gives a closer approximation. • Columns reach a rigid substratum or their length (L) is higher than critical (&2B). • Soil profile is divided into slices that are independently analysed. • Plastic strains in the soil are not considered.
• Foundation loading is under drained conditions. These simplifying assumptions impose some limitations to the analytical solution. The solution is not directly applicable to large loaded areas or stratified soils because the load distribution slope changes. Additionally, the solution is not applicable to those cases where the column does not reach a rigid substratum and it is not realistic to assume L C 2B. The solution is only valid for low to moderate foundation loads (lower than 50-100 kPa for common cases) because plastic strains in the soil are not considered.
The solution may be implemented in a spreadsheet or any programming language. A simplified calculation example and the corresponding MATLAB code are presented.
The analytical solution predicts footing settlements that agree well with numerical analyses and are usually on the safe side for most cases. Good agreement is also found for vertical displacements, vertical stresses and horizontal displacements at the soil/column interface. Additionally, the solution compares well with a smallscale laboratory test available in the literature, which confirms the validity of the solution and its hypotheses. So, the proposed solution is an analytical and accurate tool for the design of groups of stone columns beneath rigid footings.