# Estimating the Bearing Capacity of a Square Footing on a Geogrid-Stabilised Granular Layer over Clay Using the Finite Element Method

Authors: Boon Tiong Chua; Kali Prasad Nepal; Hossam Abuel-Naga
Journal: Advances in Civil Engineering
 (2023)

---

## Abstract

The bearing capacity of a square footing on an unstabilised and geogrid-stabilised granular layer over clay is estimated using finite element techniques in ABAQUS (a commercial finite element software). The model incorporates a range of soil parameters and layer thicknesses. The two-layered soil system in the model assumes a strong upper granular layer and a weak lower clay layer. The parametric results of the finite element analyses are used to develop design curves with dimensionless parameters in terms of linear gradient and the ratio of the undrained shear strength of clay to the effective vertical stress of the granular layer. The design curves are compared with the published results in the literature and found to be promising. The developed design curves are applied to estimate the minimum thickness of the granular layer and are compared with other existing methods. This application illustrates the simplicity of the proposed technique and reveals that the proposed approach yields a reasonable prediction for the bearing capacity when compared with other established methods.


---

## Body

## 1. Introduction

The use of an unbound granular fill over soft clay is a common engineering solution to overcome excessive settlement and improve bearing capacity. This conventional technique sometimes leads to very thick granular fill, which may not be economical. Other traditional treatments, such as replacing the soft soil with stronger granular fill or lime stabilisation, are typically used in soft ground conditions. However, these techniques may not always be economical and/or environmentally friendly. The use of a geogrid-stabilised granular layer can be an economical and sustainable alternative for the construction of temporary roads and working platforms on soft ground. The bearing capacities can be retained while the granular layer thickness can be reduced between 30% and 50% when the granular layer is stabilised with geogrid [1].There have been extensive research studies to predict the bearing capacity of a layered soil for nearly five decades [2–17]. The earlier bearing capacity estimation techniques were developed from both experimental investigations and analytical approaches, which led to semiempirical methods. In engineering practice, two semiempirical methods based on the equilibrium concept are widely used to calculate the bearing capacity of footings on granular layers overlying clay: (i) projected area or load spread methods and (ii) punching shear methods (refer to Chua and Nepal [18] for a detailed discussion and an examination of the limitations of these two approaches). The major shortcoming of the load spread methods is the estimation or selection of the load spread angle (α) whereas the punching shear methods have difficulty determining the angle of internal friction (ϕ′) under progressive failure conditions due to ϕ′ being stress dependent, estimating the forces along the vertical shear planes, and obtaining empirical values of the coefficient of punching shear (Ks) for a wide range of soil parameters. It is worthy to note that, for large strain during progressive failure induced by loading, the soil shear strength is not constant at all points of the shear plane for dense coarse-grained soils (i.e., sand) that are characterised by strain-softening behaviour [19].This study aims to develop a numerical modelling technique to address these shortcomings and estimate the bearing capacity of square footings on an unstabilised and geogrid-stabilised granular layer over clay. The parametric study has identified various parameters and dimensionless groups that influence the bearing capacity, which are used to develop the design curves with simple nonlinear equations. A number of assumptions are made in finite element modelling: (i) the clay is homogenous with infinite depth, (ii) the geogrid is placed horizontally at the interface of the granular and clay layers, (iii) the granular layer is homogenous withγ and ϕ′ and is compacted well to a minimum thickness of 150 mm to support the square footing before clay fails, (iv) geogrids meet the industrial requirements, for example, Austroads [20] and AASHTO [21], to comply with the geosynthetic’s survivability requirements, (v) the model is applicable for shallow footing with D < B, (vi) a settlement of 20 mm is used to generate design curves, and (vii) the influence of matric suction (capillary stresses) and unsaturated soil conditions are excluded. These assumptions are not uncommon in any bearing capacity estimation technique.


## 2. Problem Definition

The bearing capacity problem investigated in this article is illustrated in Figure1. A square footing (B/L = 1, width B and length L) rests on a granular layer of thickness H, a unit weight γ, and an internal friction angle ϕ′ and carries a surcharge q. The granular layer is overlying clay with an infinite depth of undrained shear strength cu, and the deep groundwater table is below the granular layer. In this study, the short-term stability of the footing is considered, and hence, it is reasonable to assume that the granular layer is fully drained and the clay bed is in an undrained condition. Both granular and clay layers are assumed to be homogenous.Figure 1
Schematic diagram and terminology of the problem.The soil particles move in radial directions (not in a parallel plane as in strip footing) when the soil support of the square and circular footing yields [22]. Hence, a two-dimensional (2D) axisymmetric model is appropriate to model circular and square footings. Using a square area in a 2D axisymmetric model may cause mesh generation problems, especially when the radius of the model is not constant. The square footing can be translated into circular footing of an equivalent area using the conversion factor of 4/π=1.13 or vice versa [16]. It is for this reason that circular footings are used in the finite element modelling to generate responses for the two-layer soil system. These results are then converted into square feet to develop design curves. Three-dimensional (3-D) modelling could replicate square footing more accurately, but computation time is huge as compared to a 2D axisymmetric model. Past studies [23, 24] have found a maximum difference of around 10% between the 3D model and the 2D axisymmetric model for compressive strains and stresses in subgrade. Given that long computation times are required for a 3D model with marginal benefit, a 2D axisymmetric model is used in this study.The limiting bearing capacity for a strong granular layer over a weak clay soil system becomes constant when the critical depth-to-width ratio (H/B) is around 2.5 to 2.7 [25–27]. The finite element model (FEM) solutions in this study have been computed for problems where H/B ranges from 0.5 to 2.0 and cu varies from 15 to 120 kPa with a constant footing dimension, B. The proposed model geometry and the soil strength ranges provide coverage for most practical problems.The projected area with a constant load dispersion angleα in the granular layer from the footing to the clay surface translates into a rigid stress block of truncated cone or trapezoidal shape, as shown in Figure 1. The bearing capacity of the circular footing is determined from the ultimate bearing stress of the underlying clay and the base area of the block as follows [2]:(1)qu=1+2HBtanα2cuNcSc+p0≤qg.The bearing capacity of a shallow square footing on a clay layer can be expressed as follows [3, 16]:(2)qc=cuNcSc+p0,where p0=γD+γH and Sc=1+0.2B/L=1.2The bearing capacity factor,Nc, for the underlying clay layer varies nonlinearly with mobilised shear up to a maximum of 6.2 where no outward support is required [28–31]. In the event of a high groundwater table in the granular soil layer, the effect of groundwater table in the granular layer can be accounted for by subtracting the pore water pressure from the total pressure to obtain the effective vertical stress (po) at the base of granular layer. In addition, the effective vertical stress (po) is inclusive of overburden or surcharge stress where applicable.For two-layer soil systems with a dense granular layer overlying soft clay, equation (2) can be substituted into equation (1) to form(3)qu=1+2HBtanα2qc≤qg,where qg is the ultimate bearing capacity of the granular layer of infinite depth without clay.Equation (3) infers that the increment of bearing capacity is limited to a certain value of H/B after which it becomes constant, and the overall ultimate bearing capacity of a two-layered soil system qu cannot exceed the bearing capacity provided entirely by the thick granular layer qg [7]. This implies that the bearing capacity for a limited depth of granular layer overlying a weaker clay layer does not exceed that of a deep granular bed (where H/B = ∞).In the case of no embedment,qg can be determined from [7](4)qg=0.5BγNγSγ.Divide all the items in equation (3) by qc and the equations can be written in the form(5)quqc=1+2HBtanα2≤qgqc,(6)quqc0.5=1+mHB≤qgqc,where m=2tanα is slope or gradient.Equation (5) infers that the value of bearing capacity ratio qu/qc depends on α and H/B. This dimensionless group has a significant influence on the bearing capacity of the granular layer. When qu/qg equals unity, the bearing capacity is equal to the bearing capacity of the granular layer alone, and a critical H/B is reached at which the clay has no effect on the bearing capacity. Equation (5) can be converted into a linear equation as in equation (6). It shows that qu/qc0.5 is linearly related to H/B with m being the gradient (slope) and the intercept on the qu/qc0.5 axis at 1. This is the fundamental basis for this approach to estimating the bearing capacity of square footings on granular layer overlying clay.


### 2.1. Validation of the Linear Relationship

The proposed equation (6) showing the linear relationship of qu/qc0.5 and H/B is investigated using the published results in literature. Extensive past research data from the Load Spread Methods, Laboratory Tests (LT), Limit Analysis (denoted AT), Finite Element Analysis (FEA), Finite Element Limit Analyses (FELA), and Finite Element Methods (FEM) for the bearing capacity of square and circular footings on granular layer overlying clay have been collated and processed to test the linear correlation. The qu/qc0.5 and H/B values determined from past research data for various cu values with different ϕ′ are plotted in Figures 2 to 10 [3, 4, 6, 7, 32–40]. The processed data exhibits a strong positive linear association with points that are close to the linear trend line, as explained by higher R2 values. The m value for square footing can be determined by dividing the m value for circular footing by 1.13 [16]. This conversion factor is derived directly using the same area geometry for square and circular footings for a given applied vertical load.Figure 2
qu/qc0.5 versus H/B for an unstabilised granular layer (data source: Meyerhof [3]): (a) MAT ϕ′ = 35°, (b) MAT ϕ′ = 40°, and (c) MAT ϕ′ = 45°.
(a)(b)(c)Figure 3
qu/qc0.5 versus H/B for an unstabilised granular layer (data source: Hannah and Meyerhof [4]): (a) HM ϕ′ = 40° and (b) HM ϕ′ = 45°.
(a)(b)Figure 4
qu/qc0.5 versus H/B for an unstabilised granular layer (data source: (a) Okamura et al. [7] and (b) Laman et al. [32]): (a) LT ϕ′ = 38–40° and (b) LT ϕ′ = 43°.
(a)(b)Figure 5
qu/qc0.5 versus H/B for an unstabilised granular layer (data source: (a) Lee et al. [33] and (b) Lee et al. [34]): (a) LT ϕ′ = 37°–39° and (b) AT ϕ′ = 37°–39°.
(a)(b)Figure 6
qu/qc0.5 versus H/B for an unstabilised granular layer (data source: (a) Lee et al. [33] and (b) Hemalatha et al. [35]): (a) LT ϕ′ = 31° and (b) LT ϕ′ = 43°.
(a)(b)Figure 7
qu/qc0.5 versus H/B for unstabilised granular layer (data source: Kumar and Chakraborty [36]): (a) FELA ϕ′ = 40° and (b) FELA ϕ′ = 37°–39°.
(a)(b)Figure 8
qu/qc0.5 versus H/B for an unstabilised granular layer (data source: Park and Park [37]): (a) FEM ϕ′ = 35°, (b) FEM ϕ′ = 40°, and (c) FEM ϕ′ = 45°.
(a)(b)(c)Figure 9
qu/qc0.5 versus H/B for an unstabilised granular layer (data source: (a) Gupta et al. [38] and (b) load spread method).
(a)(b)Figure 10
qu/qc0.5 versus H/B for geogrid-stabilised granular layer (data source: (a) Lees and Matthias [39] and (b) Roy and Deb [40]): (a) FEA ϕ′ = 50° and (b) LT ϕ′ = 38°.
(a)(b)By assuming the load spread angles of 2V : 1H (26.6°) and 1V : 1H (45°), the bearing capacity of square footing on a granular layer over soft clay is calculated by using the classical load spread angle method. The results are then used to generatequ/qc0.5 and H/B values. The load spread method results show a strong linear relationship between qu/qc0.5 and H/B as illustrated in Figure 9(b).In all cases,qu/qc0.5 starts at 1 for H/B = 0 when the footing is directly on the subgrade. As the H/B increases, so does the contribution of the granular layer to the bearing capacity, according to qu/qc0.5. The qu/qc0.5 value is higher for low clay strength cu because shear stresses are relatively low in the granular layer and it can maintain a higher load spread. In the case of higher clay strength cu values, the granular layer shear stresses are higher at bearing capacity failure and the load spread is narrower, resulting in a more concentrated load with depth. The lines shown in Figures 2 to 10 will eventually reach the critical depth as H/B increases when qu=qg beyond which the shear failure occurs entirely within the granular layer rather than punching shear failure becoming critical. However, none of the source data has covered beyond this test limit otherwise a minimum ultimate bearing capacity for one of the soil layers govern because the soil will fail at the lowest failure load encountered.The rigorous assessment of published data from various past studies clearly demonstrates that there is a clear indication of a linear relationship betweenqu/qc0.5 and H/B with intercept at 1 for square footing resting on unstabilised and geogrid-stabilised granular layer over clay. Hence, it is reasonable to conclude that the linear relationship is value and can be used.


### 2.2. Finite Element Model

The square shape of a loading area can be translated into a circular shape or vice versa [41]. The soil particles move in radial directions (not in a parallel plane) when the soil support of the square and circular footings yields. The circular area can be used as a representative contact area in axisymmetric modelling [24, 42–46]. A 2D axisymmetric model is developed using a commercially available finite element program, ABAQUS, as shown in Figure 11. This 2D axisymmetric model is selected because of its simplicity, greater efficiency regarding computation time execution time, ease of entry and reduction of data, and adaptability for parametric study.Figure 11
Axisymmetric FEM model geometry and boundary conditions with circular loading on a granular layer over a clay bed.The dimensions of the model geometry are selected to minimise the undesirable resonance effects of the boundary’s reflection. The geogrid reinforcement is placed between the clay and granular layers. The conventional kinematic boundary conditions are adopted when the vertical sides are fixed in a horizontal direction with vertical displacement and the bottom of the model is encastred (fixed) to both vertical and horizontal directions. The footing is assumed to be in perfect contact with the granular base (that is, they share the nodal points in between). The static load is applied using the vertical displacement method up to 20 mm above the footing widthB to simulate the mechanical stabilisation of the granular layer. This implies the retention of the granular layer thickness and geometry throughout its life. The geogrid-tensioned membrane effect is not applicable without the development of large surface ruts.Eight node biquadratic axisymmetric quadrilateral solid elements and a reduced integration technique (element type CAX8R) were implemented to discretise the subgrade and granular base, while twonode linear axisymmetric membrane elements (MAX1) were used to discretise the geogrid reinforcement with a thickness of 0.003 m [47].For the purpose of this mesh convergence study, mesh sizes with constant topography ranging from 30 mm to 100 mm were tested on the whole system, with responses to key indicators. The comparisons of the results on surface deformation, vertical strain, and vertical compressive stress at the top of the subgrade obtained for each mesh size are summarised in Table1. As the mesh size becomes finer, the system performance indicators converge into a relatively constant value. It is clear that a 30 mm mesh size will minimise any discretisation errors and achieve a good approximation to the exact solution. Although the 40 mm size mesh has less than half the number of elements compared to the fine mesh (30 mm), the results are comparatively close to those of the fine mesh. Based on the consideration of accuracy, computation time, and memory storage for the finite element simulations, the optimum equivalent element size was found to be 40 mm. However, the adopted equivalent 40 mm finite element mesh topography is refined to make it well graded, with the smaller elements concentrated under the loaded area and coarser mesh further from the loaded area, as illustrated in Figure 12(a). The selected mesh density was therefore used for the remaining analyses. It is worth noting that to reduce the numerical problems due to mesh dependency, which arise when soils with strain-softening behaviour are involved in the deformation processes, the methodology proposed by Troncone et al. [48] can be adopted.Table 1
Summary of mesh convergence results.
Equivalent mesh size (mm)Total number of elementsδTop of granular layer(mm)σvTop of clay(kPa)ɛvTop of clay(μɛ)302100013.3182.30.03240967513.6190.60.03350756014.4168.20.028100189017.3151.80.023Figure 12
(a) 2D axisymmetric finite element mesh and (b) typical interface slip.
(a)(b)The Mohr–Coulomb Models (MCM) and Tresca Models are selected to represent elastoplastic material with hardening behaviour of granular and clay layers, respectively, whereby the shear strength developed in soil is a function of internal friction’s angle, cohesion and applied stress [49–52]. In addition, the cohesion in the MCM for granular bases is set at 0.1 kPa to enhance the stability of numerical analysis [33, 37].The influence of the associated flow rule is significant for higher friction angles (ϕ′ > 30°) up to 15%, as demonstrated by past researchers [53–56]. The Mohr–Coulomb model accounts for soil dilatancy and allows the yield surface to expand based on plastic straining. In this article, the associated flow rule is replicated by selecting the dilatancy angle ψ equal to the internal friction angle less 30° [37, 57, 58]. The modulus of the granular base is dependent on the modulus of the underlying clay. Soft clay soil does not provide the support needed to obtain good compaction of the granular base. The modulus ratio of the granular base to the clay subgrade soil is usually taken to be 1.0 to 5.0 as recommended by Heukelom and Klomp [59], and a ratio of 2.0 to 5.0 was selected for the FEM analysis.For 20 mm, footing settlement, the induced strain in the geogrid is very small and is considered to be within the elastic range. Therefore, a linear elastic model is used to characterise the behaviour of the geogrid material. In addition, the maximum geogrid tensile strength (axial stiffness) at 2% strain from the material in the machine and cross-machine directions is used. Other researchers [60–64] have successfully used the same model assumption. The mechanical stabilising effect of geogrid extends over 300 mm from the geogrid plane is simulated in the model by assigning stress concentrations of 20 kPa at the geogrid elevation and decreases linearly at a constant rate to nominal 0.1 kPa at 300 mm above the geogrid plane within the influence zone [15]. The geogrid selection is based on commonly used in engineering practices. Other presumptive material constitutive models input parameters considered in this parametric study are based on published data from Ameratung et al. [58], BRE [65], and Erickson and Drescher [66] as summarised in Table 2.Table 2
Material properties for FEM analysis of square footings.
MaterialsElement (type)Model and parametersModulus (MPa)Poisson ratio (ν)Unit weight (kg/m3)CommentsGranular baseSolid (CAX8R)Mohr–Coulombϕ″ = 35°–45°c′ = 0.1 kPa10–2000.351900–2100Loose to dense granular baseψ = ϕ′–30°ClaySolid (CAX8R)Mohr–Coulombcu = 15–150 kPa5–1000.451700–2000Very soft to stiff clayGeogridMembrane (MAX1)Elastic2000.35100–340Axial stiffness 600 kN/mNote. Material parameters based on compilation of published data from Ameratung et al. [58], BRE [65], Erickson and Drescher [66] and global synthetics technical [67].The shear resistance interaction between geogrid reinforcement and the granular base is replicated using the Coulomb friction model by inputting two material properties, namely, friction coefficient and elastic slip. A small elastic slippage of 0.001 m is likely to occur for a small settlement of 20 mm [50, 61, and 68], as illustrated in Figure 12(b). The friction coefficients for geogrid-granular base and geogrid-clay are 0.84 and 0.62, respectively are based on published test results in BOSTD [69].


## 2.1. Validation of the Linear Relationship

The proposed equation (6) showing the linear relationship of qu/qc0.5 and H/B is investigated using the published results in literature. Extensive past research data from the Load Spread Methods, Laboratory Tests (LT), Limit Analysis (denoted AT), Finite Element Analysis (FEA), Finite Element Limit Analyses (FELA), and Finite Element Methods (FEM) for the bearing capacity of square and circular footings on granular layer overlying clay have been collated and processed to test the linear correlation. The qu/qc0.5 and H/B values determined from past research data for various cu values with different ϕ′ are plotted in Figures 2 to 10 [3, 4, 6, 7, 32–40]. The processed data exhibits a strong positive linear association with points that are close to the linear trend line, as explained by higher R2 values. The m value for square footing can be determined by dividing the m value for circular footing by 1.13 [16]. This conversion factor is derived directly using the same area geometry for square and circular footings for a given applied vertical load.Figure 2
qu/qc0.5 versus H/B for an unstabilised granular layer (data source: Meyerhof [3]): (a) MAT ϕ′ = 35°, (b) MAT ϕ′ = 40°, and (c) MAT ϕ′ = 45°.
(a)(b)(c)Figure 3
qu/qc0.5 versus H/B for an unstabilised granular layer (data source: Hannah and Meyerhof [4]): (a) HM ϕ′ = 40° and (b) HM ϕ′ = 45°.
(a)(b)Figure 4
qu/qc0.5 versus H/B for an unstabilised granular layer (data source: (a) Okamura et al. [7] and (b) Laman et al. [32]): (a) LT ϕ′ = 38–40° and (b) LT ϕ′ = 43°.
(a)(b)Figure 5
qu/qc0.5 versus H/B for an unstabilised granular layer (data source: (a) Lee et al. [33] and (b) Lee et al. [34]): (a) LT ϕ′ = 37°–39° and (b) AT ϕ′ = 37°–39°.
(a)(b)Figure 6
qu/qc0.5 versus H/B for an unstabilised granular layer (data source: (a) Lee et al. [33] and (b) Hemalatha et al. [35]): (a) LT ϕ′ = 31° and (b) LT ϕ′ = 43°.
(a)(b)Figure 7
qu/qc0.5 versus H/B for unstabilised granular layer (data source: Kumar and Chakraborty [36]): (a) FELA ϕ′ = 40° and (b) FELA ϕ′ = 37°–39°.
(a)(b)Figure 8
qu/qc0.5 versus H/B for an unstabilised granular layer (data source: Park and Park [37]): (a) FEM ϕ′ = 35°, (b) FEM ϕ′ = 40°, and (c) FEM ϕ′ = 45°.
(a)(b)(c)Figure 9
qu/qc0.5 versus H/B for an unstabilised granular layer (data source: (a) Gupta et al. [38] and (b) load spread method).
(a)(b)Figure 10
qu/qc0.5 versus H/B for geogrid-stabilised granular layer (data source: (a) Lees and Matthias [39] and (b) Roy and Deb [40]): (a) FEA ϕ′ = 50° and (b) LT ϕ′ = 38°.
(a)(b)By assuming the load spread angles of 2V : 1H (26.6°) and 1V : 1H (45°), the bearing capacity of square footing on a granular layer over soft clay is calculated by using the classical load spread angle method. The results are then used to generatequ/qc0.5 and H/B values. The load spread method results show a strong linear relationship between qu/qc0.5 and H/B as illustrated in Figure 9(b).In all cases,qu/qc0.5 starts at 1 for H/B = 0 when the footing is directly on the subgrade. As the H/B increases, so does the contribution of the granular layer to the bearing capacity, according to qu/qc0.5. The qu/qc0.5 value is higher for low clay strength cu because shear stresses are relatively low in the granular layer and it can maintain a higher load spread. In the case of higher clay strength cu values, the granular layer shear stresses are higher at bearing capacity failure and the load spread is narrower, resulting in a more concentrated load with depth. The lines shown in Figures 2 to 10 will eventually reach the critical depth as H/B increases when qu=qg beyond which the shear failure occurs entirely within the granular layer rather than punching shear failure becoming critical. However, none of the source data has covered beyond this test limit otherwise a minimum ultimate bearing capacity for one of the soil layers govern because the soil will fail at the lowest failure load encountered.The rigorous assessment of published data from various past studies clearly demonstrates that there is a clear indication of a linear relationship betweenqu/qc0.5 and H/B with intercept at 1 for square footing resting on unstabilised and geogrid-stabilised granular layer over clay. Hence, it is reasonable to conclude that the linear relationship is value and can be used.


## 2.2. Finite Element Model

The square shape of a loading area can be translated into a circular shape or vice versa [41]. The soil particles move in radial directions (not in a parallel plane) when the soil support of the square and circular footings yields. The circular area can be used as a representative contact area in axisymmetric modelling [24, 42–46]. A 2D axisymmetric model is developed using a commercially available finite element program, ABAQUS, as shown in Figure 11. This 2D axisymmetric model is selected because of its simplicity, greater efficiency regarding computation time execution time, ease of entry and reduction of data, and adaptability for parametric study.Figure 11
Axisymmetric FEM model geometry and boundary conditions with circular loading on a granular layer over a clay bed.The dimensions of the model geometry are selected to minimise the undesirable resonance effects of the boundary’s reflection. The geogrid reinforcement is placed between the clay and granular layers. The conventional kinematic boundary conditions are adopted when the vertical sides are fixed in a horizontal direction with vertical displacement and the bottom of the model is encastred (fixed) to both vertical and horizontal directions. The footing is assumed to be in perfect contact with the granular base (that is, they share the nodal points in between). The static load is applied using the vertical displacement method up to 20 mm above the footing widthB to simulate the mechanical stabilisation of the granular layer. This implies the retention of the granular layer thickness and geometry throughout its life. The geogrid-tensioned membrane effect is not applicable without the development of large surface ruts.Eight node biquadratic axisymmetric quadrilateral solid elements and a reduced integration technique (element type CAX8R) were implemented to discretise the subgrade and granular base, while twonode linear axisymmetric membrane elements (MAX1) were used to discretise the geogrid reinforcement with a thickness of 0.003 m [47].For the purpose of this mesh convergence study, mesh sizes with constant topography ranging from 30 mm to 100 mm were tested on the whole system, with responses to key indicators. The comparisons of the results on surface deformation, vertical strain, and vertical compressive stress at the top of the subgrade obtained for each mesh size are summarised in Table1. As the mesh size becomes finer, the system performance indicators converge into a relatively constant value. It is clear that a 30 mm mesh size will minimise any discretisation errors and achieve a good approximation to the exact solution. Although the 40 mm size mesh has less than half the number of elements compared to the fine mesh (30 mm), the results are comparatively close to those of the fine mesh. Based on the consideration of accuracy, computation time, and memory storage for the finite element simulations, the optimum equivalent element size was found to be 40 mm. However, the adopted equivalent 40 mm finite element mesh topography is refined to make it well graded, with the smaller elements concentrated under the loaded area and coarser mesh further from the loaded area, as illustrated in Figure 12(a). The selected mesh density was therefore used for the remaining analyses. It is worth noting that to reduce the numerical problems due to mesh dependency, which arise when soils with strain-softening behaviour are involved in the deformation processes, the methodology proposed by Troncone et al. [48] can be adopted.Table 1
Summary of mesh convergence results.
Equivalent mesh size (mm)Total number of elementsδTop of granular layer(mm)σvTop of clay(kPa)ɛvTop of clay(μɛ)302100013.3182.30.03240967513.6190.60.03350756014.4168.20.028100189017.3151.80.023Figure 12
(a) 2D axisymmetric finite element mesh and (b) typical interface slip.
(a)(b)The Mohr–Coulomb Models (MCM) and Tresca Models are selected to represent elastoplastic material with hardening behaviour of granular and clay layers, respectively, whereby the shear strength developed in soil is a function of internal friction’s angle, cohesion and applied stress [49–52]. In addition, the cohesion in the MCM for granular bases is set at 0.1 kPa to enhance the stability of numerical analysis [33, 37].The influence of the associated flow rule is significant for higher friction angles (ϕ′ > 30°) up to 15%, as demonstrated by past researchers [53–56]. The Mohr–Coulomb model accounts for soil dilatancy and allows the yield surface to expand based on plastic straining. In this article, the associated flow rule is replicated by selecting the dilatancy angle ψ equal to the internal friction angle less 30° [37, 57, 58]. The modulus of the granular base is dependent on the modulus of the underlying clay. Soft clay soil does not provide the support needed to obtain good compaction of the granular base. The modulus ratio of the granular base to the clay subgrade soil is usually taken to be 1.0 to 5.0 as recommended by Heukelom and Klomp [59], and a ratio of 2.0 to 5.0 was selected for the FEM analysis.For 20 mm, footing settlement, the induced strain in the geogrid is very small and is considered to be within the elastic range. Therefore, a linear elastic model is used to characterise the behaviour of the geogrid material. In addition, the maximum geogrid tensile strength (axial stiffness) at 2% strain from the material in the machine and cross-machine directions is used. Other researchers [60–64] have successfully used the same model assumption. The mechanical stabilising effect of geogrid extends over 300 mm from the geogrid plane is simulated in the model by assigning stress concentrations of 20 kPa at the geogrid elevation and decreases linearly at a constant rate to nominal 0.1 kPa at 300 mm above the geogrid plane within the influence zone [15]. The geogrid selection is based on commonly used in engineering practices. Other presumptive material constitutive models input parameters considered in this parametric study are based on published data from Ameratung et al. [58], BRE [65], and Erickson and Drescher [66] as summarised in Table 2.Table 2
Material properties for FEM analysis of square footings.
MaterialsElement (type)Model and parametersModulus (MPa)Poisson ratio (ν)Unit weight (kg/m3)CommentsGranular baseSolid (CAX8R)Mohr–Coulombϕ″ = 35°–45°c′ = 0.1 kPa10–2000.351900–2100Loose to dense granular baseψ = ϕ′–30°ClaySolid (CAX8R)Mohr–Coulombcu = 15–150 kPa5–1000.451700–2000Very soft to stiff clayGeogridMembrane (MAX1)Elastic2000.35100–340Axial stiffness 600 kN/mNote. Material parameters based on compilation of published data from Ameratung et al. [58], BRE [65], Erickson and Drescher [66] and global synthetics technical [67].The shear resistance interaction between geogrid reinforcement and the granular base is replicated using the Coulomb friction model by inputting two material properties, namely, friction coefficient and elastic slip. A small elastic slippage of 0.001 m is likely to occur for a small settlement of 20 mm [50, 61, and 68], as illustrated in Figure 12(b). The friction coefficients for geogrid-granular base and geogrid-clay are 0.84 and 0.62, respectively are based on published test results in BOSTD [69].


## 3. Results and Discussion

### 3.1. Verification of the Finite Element Model

The developed axisymmetric finite element model is validated by comparing the results with published data reported by Leng and Gabr [70]. All the input data comprising the loading conditions, constitutive models, geogrid reinforcements, and friction coefficient properties at granular base-geogrid-subgrade interfaces adopted are similar. This thereby eliminates any variations and assumptions which may have inherent inaccuracies. The outcome of this validation process provides a level of confidence in the FEM.By varying the granular layer thicknessH, the FEM is assessed for the critical responses covering the surface deformation δ of the granular layer beneath the centre of the circular footing and the vertical compressive stress σ on top of the clay. The FEM results of the critical responses are generally in good agreement, and the differences are systematic with the published data from Leng and Gabr [70] for both the unstabilised and geogrid-stabilised granular layer over clay, as demonstrated in Figure 13.Figure 13
Surface deformations for unstabilised and geogrid-stabilised pavement structures (US: unstabilised and GS: geogrid-stabilised).
(a)(b)


### 3.2. Proposed Design Curves for Square Footings on Granular Layer over Clay

Them values derived from the linear relationship of (qu/qc)0.5 and H/B for unstabilised and geogrid-stabilised granular layers have been used to formulate the design curves. A lower bound fitted with power trend line option is used to express the relationship of m and cu/p0 to accommodate the variability of data. This approach demonstrates the simple derivation of the bearing capacity directly from the shear strengths of the individual soil layers and makes it easy to visualise the contribution of the individual strengths to the overall bearing capacity, as shown in Figure 14.Figure 14
Design curves for square footing on an unstabilised and geogrid-stabilised granular layer over clay. (a) Nonstabilised. (b) Geogrid-stabilised.
(a)(b)The nonlinear expression for unstabilised and geogrid-stabilised granular base with various values of internal friction angleϕ′ ranging from 35° to 45° is shown in Figure 14. The design curves with nonlinear expression allow a direct determination of m for a given ϕ′ and γ of the granular layer and cu/p0 at the base of granular layer graphically or by using the appropriate nonlinear equations as shown in Figure 14. The m value is then applied to (6) to determine the bearing capacity of the square footing on an unstabilised or geogrid-stabilised granular layer over clay. For example, the value of m for an unstabilised granular layer (H = 0.37 m, ϕ′ = 40° and γ = 20 kN/m3) over clay (cu = 40 kN/m2) from Figure 14(a) is 0.34, and qc in (6) is equal to 1.2 × 6.2 × 40 = 297.6 kPa for a square footing of 0.5 m by 0.5 m. Substituting these values in (6) results bearing capacity, qu = 468.5 kPa.The design curves for unstabilised and geogrid-stabilised granular bases show a similar trend, with a significant shift towards a higherm value for all cases of cu/p0 as illustrated in Figure 14. The m value is very sensitive when the shear strength of clay is very low. This is evident with the exponential increase in the m value when cu/p0 is less than approximately 2.5. This can be attributed to the unstable nature of very soft clay.


### 3.3. Comparison with Published Data

The validity of the proposed method for the estimation of bearing capacity is examined by verifying the applicability of this technique in comparison with other established methods data. Whenm values derived from the linear relationship of (qu/qc)0.5 and H/B for an unstabilised granular layer using a variety of methods are superimposed on the design curves, it shows that there are some variabilities, suggesting some dependency on other variables. Therefore, it is reasonable to adopt lower bound power trend lines for various ϕ′ values for the intended practical design use as depicted in Figure 15.Figure 15
m values versus cu/po from various studies for an unstabilised granular layer with ϕ′ = 35°, 40°, and 45°.The design chart is compared with the load spread method for square footings with assumed load spread angles of 26.6° (2V : 1H) and 45° (1V : 1H) and the correspondingm values of 0.9 and 1.8, as shown in Figure 16(a). It can be seen that the m values derived from the load spread method always remain on the higher side in comparison to the results of other methods. In comparison with other available models, a smaller load spread angle will yield comparable results. At lower cu/p0 values, the m values show an increase because of the increased influence of granular layer self-weight at lower cu values in equation (1). This comparison clearly reveals that the adoption of a single angle of load spread in all cases could result in over or underestimation of bearing capacity.Figure 16
m values versus cu/po from various studies for (a) unstabilised and (b) geogrid-stabilised granular layer with ϕ′ = 35°, 40°, 45°, and 50°.
(a)(b)In the case of a geogrid-stabilised granular layer, the design curves indicate a similar trend when compared with Lee’s [15] published data. Figure 16(b) shows that the nonlinear relationships proposed by Lees and Matthias [39] are shown to have higher m values for all cases of cu/p0 than the design curves because of the denser geogrid-stabilised granular layer with an internal friction angle ϕ′ = 50° as compared with ϕ′ = 35° to 45° being used in this study.


### 3.4. Comparative Evaluation of Minimum Granular Thickness

The use of the developedm and cu/p0 design curves for square footing loaded on the surface of unstabilised and geogrid-stabilised granular layers without surface surcharge to calculate the minimum thickness of a granular layer is illustrated by a typical case. The purpose of this case is to show the application of the design curves in the design of square footing and to compare against established design approaches, such as Meyerhof [3], Hanna and Meyerhof [4], load spread, Lees [15, 16], and Lees and Matthias [39] methods to estimate the minimum thickness of an unstabilised and geogrid-stabilised granular layer. The input soil parameters and loading conditions for the case are summarised in Table 3 and the results are plotted in Figure 17.Table 3
Soil parameters and loading conditions.
MaterialsSoil parametersUnit weight (kN/m3)Loading conditionsGranular fillϕ″ = 40°, c′ = 0 kPa20Surface loading = 400 kN/m2Claycu = 20 kPa19Footing dimensions: 0.5 m by 0.5 mFigure 17
Granular layer thickness versus clay undrained shear strength (US: unstabilised and GS: geogrid-stabilised).As shown in Figure17, by using α equal to 26.6° (2V : 1H) and 45° (1V : 1H), the estimated bearing capacities are higher than other methods. The granular layer thickness decreases with the increase in α. This can be attributed to the decrease in load transfer to the interface resulting in the increase in ultimate bearing capacity, and vice versa. This shows that the α has a significant influence on the vertical compressive stress, which in turn affects the bearing capacity. This observation explains the rationale for using a smaller α between 11.3° (5V : 1H) to 18.4° (3V : 1H) as recommended in practice [12]. The results in this study for the unstabilised granular base have good agreement with the T-value method by Lees [16]. The bearing capacities are higher when compared with the established classical methods by Meyerhof [3] and Hanna and Meyerhof [4] resulting in thinner unstabilised granular layer thickness. It is worth noting that the punching shear model proposed by Meyerhof [3] and Hanna and Meyerhof [4] assumes the bearing capacity of the clay is unaffected by the vertical stresses applied to the clay surface immediately outside of the loaded area. Madhav and Sharma [71] found that the surcharge load can increase the bearing capacity up to 30%. This may have been attributed to the underestimation of bearing capacity, resulting in a thicker granular layer in comparison with other methods.The load spread method assumes the vertical stresses generated from the footing is dispersed in the granular layer to the clay surface at a constant angleα to the vertical and disregards the mobilised shear resistance in granular base. It also assumes that the stress increase under the centre of the loaded area equals the stress increase under the corner or side of the loaded area. In reality, the soil underneath the centre of the loaded area is subjected to a higher vertical stress than the soil underneath the corner or edge of the loaded area. These are attributed to an overestimation of bearing capacity resulting in thinner granular base thickness is required to withstand the footing load as evident in Figure 17. In comparison with other available models, a smaller load spread angle will yield comparable results. The value of α is dependent on the shear strength of the clay, the stiffness ratio between the layers, and the model geometry relating to the depth-to-width ratio (H/B) [9, 11, 16, 72, and 73].Without the geogrid, there is no confinement effect, and the particle interlock within the granular layer is overcome more easily, resulting in a steeper punching shear plane and a smaller load to reach failure. Overall, all the methods used in this example exhibit a similar trend with the granular layer becoming thinner with increasing values ofcu.The bearing capacity of the geogrid-stabilised granular layer in this study is in good agreement with Lees and Matthias [39] T-value method for clay undrained shear strength values between 20 kPa and 40 kPa. The geogrid confines the soil particles around it, which prevent rearrangement, thereby maintaining particle interlock through confinement effect and resisting the punching shear. This results in a granular layer that distributes the footing load on the clay foundation across a wider punching shear plane and will require a higher load to reach failure. This explains why the geogrid-stabilised granular layer has a higher bearing capacity, resulting in a thinner granular layer thickness H as illustrated in this design example.At highercu values, the geogrid-stabilised granular layer used in this study may have predicted a higher bearing capacity, which translates into a thinner granular layer for this design example. The comparisons with the load spreading method show that when the clay undrained shear strength cu increases, the values of granular layer thickness H decrease in all cases and converge to a nominal or zero value. In all cases, the geogrid-stabilised granular layer is thinner than the unstabilised granular layer, except for the load spread method. This is because of aggregate confinement provided by geogrid reinforcements through interlocking and stiffens the aggregate layer resulting in overall increase in bearing capacity of the geogrid-stabilised granular layer. It should be noted that the confinement effect constrained the lateral movement of aggregate within 300 mm above the geogrid plane in the granular layer and down within the clay layer.


## 3.1. Verification of the Finite Element Model

The developed axisymmetric finite element model is validated by comparing the results with published data reported by Leng and Gabr [70]. All the input data comprising the loading conditions, constitutive models, geogrid reinforcements, and friction coefficient properties at granular base-geogrid-subgrade interfaces adopted are similar. This thereby eliminates any variations and assumptions which may have inherent inaccuracies. The outcome of this validation process provides a level of confidence in the FEM.By varying the granular layer thicknessH, the FEM is assessed for the critical responses covering the surface deformation δ of the granular layer beneath the centre of the circular footing and the vertical compressive stress σ on top of the clay. The FEM results of the critical responses are generally in good agreement, and the differences are systematic with the published data from Leng and Gabr [70] for both the unstabilised and geogrid-stabilised granular layer over clay, as demonstrated in Figure 13.Figure 13
Surface deformations for unstabilised and geogrid-stabilised pavement structures (US: unstabilised and GS: geogrid-stabilised).
(a)(b)


## 3.2. Proposed Design Curves for Square Footings on Granular Layer over Clay

Them values derived from the linear relationship of (qu/qc)0.5 and H/B for unstabilised and geogrid-stabilised granular layers have been used to formulate the design curves. A lower bound fitted with power trend line option is used to express the relationship of m and cu/p0 to accommodate the variability of data. This approach demonstrates the simple derivation of the bearing capacity directly from the shear strengths of the individual soil layers and makes it easy to visualise the contribution of the individual strengths to the overall bearing capacity, as shown in Figure 14.Figure 14
Design curves for square footing on an unstabilised and geogrid-stabilised granular layer over clay. (a) Nonstabilised. (b) Geogrid-stabilised.
(a)(b)The nonlinear expression for unstabilised and geogrid-stabilised granular base with various values of internal friction angleϕ′ ranging from 35° to 45° is shown in Figure 14. The design curves with nonlinear expression allow a direct determination of m for a given ϕ′ and γ of the granular layer and cu/p0 at the base of granular layer graphically or by using the appropriate nonlinear equations as shown in Figure 14. The m value is then applied to (6) to determine the bearing capacity of the square footing on an unstabilised or geogrid-stabilised granular layer over clay. For example, the value of m for an unstabilised granular layer (H = 0.37 m, ϕ′ = 40° and γ = 20 kN/m3) over clay (cu = 40 kN/m2) from Figure 14(a) is 0.34, and qc in (6) is equal to 1.2 × 6.2 × 40 = 297.6 kPa for a square footing of 0.5 m by 0.5 m. Substituting these values in (6) results bearing capacity, qu = 468.5 kPa.The design curves for unstabilised and geogrid-stabilised granular bases show a similar trend, with a significant shift towards a higherm value for all cases of cu/p0 as illustrated in Figure 14. The m value is very sensitive when the shear strength of clay is very low. This is evident with the exponential increase in the m value when cu/p0 is less than approximately 2.5. This can be attributed to the unstable nature of very soft clay.


## 3.3. Comparison with Published Data

The validity of the proposed method for the estimation of bearing capacity is examined by verifying the applicability of this technique in comparison with other established methods data. Whenm values derived from the linear relationship of (qu/qc)0.5 and H/B for an unstabilised granular layer using a variety of methods are superimposed on the design curves, it shows that there are some variabilities, suggesting some dependency on other variables. Therefore, it is reasonable to adopt lower bound power trend lines for various ϕ′ values for the intended practical design use as depicted in Figure 15.Figure 15
m values versus cu/po from various studies for an unstabilised granular layer with ϕ′ = 35°, 40°, and 45°.The design chart is compared with the load spread method for square footings with assumed load spread angles of 26.6° (2V : 1H) and 45° (1V : 1H) and the correspondingm values of 0.9 and 1.8, as shown in Figure 16(a). It can be seen that the m values derived from the load spread method always remain on the higher side in comparison to the results of other methods. In comparison with other available models, a smaller load spread angle will yield comparable results. At lower cu/p0 values, the m values show an increase because of the increased influence of granular layer self-weight at lower cu values in equation (1). This comparison clearly reveals that the adoption of a single angle of load spread in all cases could result in over or underestimation of bearing capacity.Figure 16
m values versus cu/po from various studies for (a) unstabilised and (b) geogrid-stabilised granular layer with ϕ′ = 35°, 40°, 45°, and 50°.
(a)(b)In the case of a geogrid-stabilised granular layer, the design curves indicate a similar trend when compared with Lee’s [15] published data. Figure 16(b) shows that the nonlinear relationships proposed by Lees and Matthias [39] are shown to have higher m values for all cases of cu/p0 than the design curves because of the denser geogrid-stabilised granular layer with an internal friction angle ϕ′ = 50° as compared with ϕ′ = 35° to 45° being used in this study.


## 3.4. Comparative Evaluation of Minimum Granular Thickness

The use of the developedm and cu/p0 design curves for square footing loaded on the surface of unstabilised and geogrid-stabilised granular layers without surface surcharge to calculate the minimum thickness of a granular layer is illustrated by a typical case. The purpose of this case is to show the application of the design curves in the design of square footing and to compare against established design approaches, such as Meyerhof [3], Hanna and Meyerhof [4], load spread, Lees [15, 16], and Lees and Matthias [39] methods to estimate the minimum thickness of an unstabilised and geogrid-stabilised granular layer. The input soil parameters and loading conditions for the case are summarised in Table 3 and the results are plotted in Figure 17.Table 3
Soil parameters and loading conditions.
MaterialsSoil parametersUnit weight (kN/m3)Loading conditionsGranular fillϕ″ = 40°, c′ = 0 kPa20Surface loading = 400 kN/m2Claycu = 20 kPa19Footing dimensions: 0.5 m by 0.5 mFigure 17
Granular layer thickness versus clay undrained shear strength (US: unstabilised and GS: geogrid-stabilised).As shown in Figure17, by using α equal to 26.6° (2V : 1H) and 45° (1V : 1H), the estimated bearing capacities are higher than other methods. The granular layer thickness decreases with the increase in α. This can be attributed to the decrease in load transfer to the interface resulting in the increase in ultimate bearing capacity, and vice versa. This shows that the α has a significant influence on the vertical compressive stress, which in turn affects the bearing capacity. This observation explains the rationale for using a smaller α between 11.3° (5V : 1H) to 18.4° (3V : 1H) as recommended in practice [12]. The results in this study for the unstabilised granular base have good agreement with the T-value method by Lees [16]. The bearing capacities are higher when compared with the established classical methods by Meyerhof [3] and Hanna and Meyerhof [4] resulting in thinner unstabilised granular layer thickness. It is worth noting that the punching shear model proposed by Meyerhof [3] and Hanna and Meyerhof [4] assumes the bearing capacity of the clay is unaffected by the vertical stresses applied to the clay surface immediately outside of the loaded area. Madhav and Sharma [71] found that the surcharge load can increase the bearing capacity up to 30%. This may have been attributed to the underestimation of bearing capacity, resulting in a thicker granular layer in comparison with other methods.The load spread method assumes the vertical stresses generated from the footing is dispersed in the granular layer to the clay surface at a constant angleα to the vertical and disregards the mobilised shear resistance in granular base. It also assumes that the stress increase under the centre of the loaded area equals the stress increase under the corner or side of the loaded area. In reality, the soil underneath the centre of the loaded area is subjected to a higher vertical stress than the soil underneath the corner or edge of the loaded area. These are attributed to an overestimation of bearing capacity resulting in thinner granular base thickness is required to withstand the footing load as evident in Figure 17. In comparison with other available models, a smaller load spread angle will yield comparable results. The value of α is dependent on the shear strength of the clay, the stiffness ratio between the layers, and the model geometry relating to the depth-to-width ratio (H/B) [9, 11, 16, 72, and 73].Without the geogrid, there is no confinement effect, and the particle interlock within the granular layer is overcome more easily, resulting in a steeper punching shear plane and a smaller load to reach failure. Overall, all the methods used in this example exhibit a similar trend with the granular layer becoming thinner with increasing values ofcu.The bearing capacity of the geogrid-stabilised granular layer in this study is in good agreement with Lees and Matthias [39] T-value method for clay undrained shear strength values between 20 kPa and 40 kPa. The geogrid confines the soil particles around it, which prevent rearrangement, thereby maintaining particle interlock through confinement effect and resisting the punching shear. This results in a granular layer that distributes the footing load on the clay foundation across a wider punching shear plane and will require a higher load to reach failure. This explains why the geogrid-stabilised granular layer has a higher bearing capacity, resulting in a thinner granular layer thickness H as illustrated in this design example.At highercu values, the geogrid-stabilised granular layer used in this study may have predicted a higher bearing capacity, which translates into a thinner granular layer for this design example. The comparisons with the load spreading method show that when the clay undrained shear strength cu increases, the values of granular layer thickness H decrease in all cases and converge to a nominal or zero value. In all cases, the geogrid-stabilised granular layer is thinner than the unstabilised granular layer, except for the load spread method. This is because of aggregate confinement provided by geogrid reinforcements through interlocking and stiffens the aggregate layer resulting in overall increase in bearing capacity of the geogrid-stabilised granular layer. It should be noted that the confinement effect constrained the lateral movement of aggregate within 300 mm above the geogrid plane in the granular layer and down within the clay layer.


## 4. Conclusions

The bearing capacity of a square footing placed on the surface of an unstabilised and geogrid-stabilised granular layer over clay has been replicated using an ABAQUS 2D axisymmetric finite element model. Extensive investigative work regarding the linearity betweenqu/qc0.5 and H/B with intercept at 1 is conducted with the help of past research studies and found to be the case. The results of the FEM analysis are found to compare reasonably well with the available published data. The finite element model is validated using past research data and found to be in good agreement. This provides confidence that the finite element modelling techniques are correct prior to conducting a parametric study to generate results for preparing the design curves.The proposed approach for the unstabilised granular layer shows good agreement when compared with theT-value method by Lees [16] and is comparable with classical semiempirical methods [3, 4] for a two-layered soil system. For the case of a geogrid-stabilised granular layer, the bearing capacity results are in good agreement with the T-value method by Lees and Matthias [39] for clay undrained shear strength values between 20 kPa and 40 kPa. Overall, the geogrid-stabilised granular layer is thinner than the unstabilised granular layer in all cases except for the load spread method.A series of dimensionless design chartsm versus cu/p0, for a wide range of cu have been developed from the modelling results. The relationship of m and cu/p0 for each granular layer’s internal friction angle ϕ′ can be expressed as nonlinear equations, which are simple and easy to use as an alternate method to design curves. The nonlinear equations can be implemented numerically using a spreadsheet. The proposed method can be applied to a square footing under embedded and surface loading similar to shallow footing (D < B) conditions, including the effect of groundwater, where the strengths of the granular layer and clay are characterised by both internal friction and cohesion. The proposed design curves or the nonlinear equations can be used for preliminary concept design to estimate granular layer thickness, as the approach is derived from a relatively small database, and it would be beneficial to verify them with field tests.



---

*Source: 3560697-2023-03-06.xml*