Stochastic numerical approach to the bearing capacity of adjacent footings on spatially variable soils

In this study, three soil models were considered to represent a range of geotechnical conditions, including both idealised and natural materials. The first two approaches represented the behaviour of natural soils but were considered in this study as theoretical assumptions. Approach 1 – perfectly cohesive soil – represented a saturated undrained clay, where only total cohesion was considered (mean value: 36.3 kPa), and the internal friction angle was assumed to be zero. Approach 2 – purely frictional soil – simulated a clean, well-compacted sand under drained conditions, with an effective friction angle (mean: 35°) and negligible cohesion (1 kPa). The cohesive-frictional soil model (Approach 3) was developed using Taranto Blue Clay (TBC), a natural marine clay deposit in southern Italy. The effective strength parameters for TBC were derived from in situ and laboratory investigations reported by Cafaro and Cherubini [14]. These include a lognormally distributed cohesion with a mean value of 36.3 kPa and a bounded friction angle with a mean of 20°, ranging from 5° to 35°. This model reflects realistic variability observed in natural cohesive-frictional soils.

Initially, the analytical calculations for bearing capacities of isolated footings were performed using Terzaghi’s formula. Two foundation widths, B = 1 and 2 m, were adopted for the calculation, and are shown in Table 1. In principle, the width B = 2 m corresponds to two foundations joined together.

Table 1

	Analytical formula	BC value for B = 1 m	BC value for B = 2 m
Approach 1 – perfectly cohesive soil	q a = 5.12 c u {q}_{a}=5.12\hspace{0.25em}{c}_{u}	185.86 kPa	185.86 kPa
Approach 2 – purely frictional soil	q a = q N q + 0.5 B γ N γ {q}_{a}=q{N}_{q}+\hspace{0.25em}0.5B\gamma {N}_{\gamma }	1108.41 kPa	1538.08 kPa
Approach 3 – cohesive-frictional soil	q a = c N c + q N q + 0.5 B γ N γ {q}_{a}=c{N}_{c}+q{N}_{q}+0.5B\gamma {N}_{\gamma }	697.43 kPa	734.77 kPa

Source: Author’s contribution.

The analytical resistances ( q a {q}_{a} ) are estimated using resistance coefficients N i {N}_{i} , according to Eurocode 7. Total cohesion was used for the cohesive soil to reflect undrained behaviour (Approach 1). In contrast, effective parameters were applied in the frictional (Approach 2) and cohesive-frictional (Approach 3) cases to represent drained conditions. This classification enables a comparative analysis of foundation performance across different assumptions about soil type. The embedment of foundations was represented by a constant load next to the foundation ( q q ) equal to 19.0 kPa 19.0\hspace{0.25em}{\rm{kPa}}\hspace{0.25em} . In the case of a perfectly cohesive soil, the load q = 0. The results of the analytical calculations, assuming mean parameter values, are presented in Table 1. The BC of purely frictional soil shows the most significant dependence on foundation width.

Prior to the stochastic simulations, a series of preliminary finite element analyses were conducted to optimise both the mesh resolution and the displacement step parameters. The domain has been discretised into 160 × 40 elements, each measuring 0.25 m × 0.25 m. The final domain dimensions were 40 m × 10 m, as illustrated in Figure 2. Such a domain of discretisation was assumed (as the effect of several tries) to ensure that no boundary effects were present. Two foundations with a dimension of B = 1 m each were adopted for the study. As in the analytical examples, the 1 m embedment was modelled by applying a pressure of 19 kN/m² to the soil surface (except in Case 1). The Mohr-Coulomb soil model was assumed for all calculations. The soil strength parameters were considered to be constant at their means presented in the previous subsection.

Figure 2

The adopted numerical model.

The elastic parameters were defined as follows: Young’s modulus E = 36,000 MPa and Poisson’s ratio υ = 0.29. This configuration enabled accurate modelling of stress redistribution between adjacent footings without artificial constraints, as confirmed by the displacement field presented in Figure 3. The figure illustrates node displacements for a deterministic case with footing spacings ranging from 0.9 to 8 m. Interaction between footings was observed for spacings up to 3 m, whereas for spacings greater than 4 m the foundations behaved independently. As evidenced by the deformation zones, the domain was sufficiently large in this case.

Figure 3

Node displacement results for a deterministic task with foundation spacing 0.9–8 m.

Additionally, a detailed sensitivity study was conducted to calibrate the displacement increment and convergence criteria for each numerical situation. The selected displacement increment step size ensured that, for large footing spacings (GAP = 8 m; Figure 3), the numerical results closely matched the analytical BC of a single footing with a width of B = 1 m (as presented in Tables 1–3). The limit state was consistently reached within 10 computational steps, thereby balancing numerical stability and physical accuracy. Numerical simulations were performed using the FEM code developed by Fenton and Griffiths [33].

Table 2

	Mean value	Standard deviation	Distribution
Cohesion (c)	36.3 kPa	19.6 kPa	Lognormal
Friction angle ( ϕ \phi )	—	—	Deterministic
Unit weight ( γ \gamma )	—	—	Deterministic

Source: Author’s contribution.

Table 3

	Mean value	Standard deviation	Distribution	Distribution parameters
Cohesion (c)	1 kPa	—	Deterministic
Friction angle ( ϕ \phi )	35°	—	Bounded	ϕ min = 25 ° {\phi }_{\min }\hspace{1em}=25^\circ
Unit weight ( γ \gamma )	19 kN/m3	—	Deterministic	ϕ max = 45 ° {\phi }_{\max }\hspace{1em}=45^\circ

Source: Author’s contribution.

The final simulations were carried out over a wide range of footing spacings denoted by GAP parameter, from 0.25, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 7.0, to 8.0 m. This range enabled a comprehensive assessment of interaction effects between adjacent foundations under various spatial configurations. The results and conclusions from this task are presented in the following sections.

All stochastic analyses were performed to describe soil strength parameters, namely, cohesion and friction angle, and following predefined statistical distributions using random fields. Random fields used in this study were generated using the Local Average Subdivision method, initially proposed by Fenton and Vanmarcke [34]. Numerical probabilistic simulations were conducted using the RFEM code developed by Fenton and Griffiths [33] and modified by the author to accommodate anisotropic correlation structures and site-specific soil variability [35]. The vertical SOF (SOF _y ) was fixed at 0.7 m. In contrast, the horizontal scale (SOF _x ) was a subject of parametric study, with values of 1, 10, and 100 m considered in the analysis. These values reflect typical depositional patterns and the geological processes that influence soil structure [36].

The selection of distributions was based on both empirical data and recent literature. Lognormal distributions were adopted for cohesion, reflecting its positively skewed nature in natural soils, as reported by field studies such as Cafaro and Cherubini [14] and by recent probabilistic analyses by Cami et al. [36]. For the friction angle, a bounded uniform distribution was used, consistent with observed variability ranges in natural deposits and contemporary stochastic modelling practices. The adopted distributions and statistical parameters for the three considered cases are summarised in Tables 2–4 for the perfectly cohesive, purely frictional, and cohesive-frictional soil models, respectively.

Table 4

	Mean value	Standard deviation	Distribution	Distribution parameters
Cohesion (c)	36.3 kPa	19.6 kPa	Lognormal
Friction angle ( ϕ \phi )	20°		Bounded	ϕ min = 5 ° {\phi }_{\min }=5^\circ
Unit weight ( γ (\gamma )	19 kN/m3	—	Deterministic	ϕ max = 35 ° {\phi }_{\max }=35^\circ

Source: Author’s contribution.

To represent spatial correlation, anisotropic Markov-type autocorrelation functions were applied, defined as: (3) ρ ( τ ) = exp − 2 τ x SO F x 2 + 2 τ y SO F y 2 . \rho (\tau )=\exp \left\{-\sqrt{{\left(\frac{2{\tau }_{x}}{{\rm{SO}}{{\rm{F}}}_{x}}\right)}^{2}+{\left(\frac{2{\tau }_{y}}{{\rm{SO}}{{\rm{F}}}_{y}}\right)}^{2}}\right\}.

Each probabilistic simulation scenario was performed using 500 realisations, a number determined through convergence testing to ensure stability of the mean BC across trials. This framework enabled a robust assessment of the influence of spatial variability on foundation performance under varying soil conditions and footing configurations.

The numerical results obtained in this study also allowed for a simple reliability analysis of adjacent footings. The results are interpreted in terms of the reliability index ( β \beta ) and compared against analytical benchmarks, providing insight into the adequacy of traditional design assumptions under realistic variability conditions. Generally, the probability of failure ( p f ) {(p}_{{\rm{f}}}) is defined as (4) p f = P [ R < S ] = Φ ( − β ) , {p}_{{\rm{f}}}=\text{P}{[}R\lt S]=\Phi (-\beta ), where R is the resistance (e.g. BC) and S is the load effect. The reliability index (β) is related to p f {p}_{{\rm{f}}} via the standard normal distribution with Φ denoting the cumulative distribution function of the standard normal distribution. Eurocode 7 does not prescribe a single target value for β but recommends values depending on the consequence class and design situation. For ultimate limit states, the reliability index corresponds to the probability of failure and is approximately p f ∼ 10 − 4 {p}_{{\rm{f}}}\sim {10}^{-4} . In present study probability of failure is assumed as (5) p f = Φ ( − β ) = P [ q < q a , d ] , {p}_{{\rm{f}}}=\Phi (-\beta )\hspace{0.25em}=\text{}P\hspace{0.25em}{[}q\lt {q}_{a,d}], where q a , d {q}_{a,d} represents the analytical BC ( q a {q}_{a} , taken from Table 1 as the reference value), divided by the safety factor (FS). Result of stochastic analysis are also presented in Section 3.

The results of the performed deterministic numerical analysis are presented in Figure 4. Across soil types, the results show pronounced differences in both the magnitude and the nature of BC responses. For perfectly cohesive soils (Figure 4a), the analytical BC remains essentially constant regardless of footing width, which is consistent with classical theory. However, numerical results indicate a slight increase as adjacent footings are brought closer together – an effect attributable to localised interaction not captured by analytical solutions. This subtle increase becomes noticeable when the distance between foundations is less than 1 m, suggesting that, even in perfectly cohesive conditions, spatial proximity can modestly enhance load BC. Still, the results are within 1% of the analytical solution, which raises doubts about whether the differences are not solely due to numerical error.

Figure 4

Analytical and numerical results for deterministic values of soil parameters for each approach: (a) ideal cohesive, (b) purely frictional, and (c) TBC.

In purely frictional soils (Figure 4b), the analytical model predicts a strong dependence on footing width, with BC increasing significantly for wider foundations (as shown in Table 1). Nevertheless, numerical simulations show that the BC of two closely spaced footings does not exceed that of a single, wider footing, which may potentially arise from partial overlap of stress zones beneath adjacent foundations. The most substantial impact of foundation width on BC, and the most significant increase due to proximity, is observed when foundations are located closer than 1.5 m to each other. Notably, even at a minimal spacing of 0.25 m, the BC does not match that of a single foundation with B = 2 m, underscoring the limitations of simple superposition in frictional soils. This interaction effect, in which proximity increases BC, is consistent with the findings of Stuart [2] and Rybak and Konderla [11] and supports the validity of the numerical model used in the present study.

For the natural cohesive-frictional TBC (Figure 4c), the behaviour is more complex. Both the width of the foundation and the GAP between adjacent foundations significantly influence BC. Remarkably, when foundations are placed less than 1 m apart, the average BC can exceed that of a single foundation with double the width (B = 2 m). Numerical results for closely spaced footings can even exceed the analytical value for a single isolated footing, suggesting a synergistic interaction effect that is entirely neglected in deterministic analytical models.

Differences preclude a direct comparison of absolute bearing capacities across soil types and foundation geometries. Figure 5a–c presents the normalised charts from Figure 4 to facilitate comparison of results across soil types. By expressing all results relative to the analytical BC of a single footing of width B = 1 m for each soil type, the normalised plots enable a direct, dimensionless comparison of interaction effects across fundamentally different soils. This approach isolates the influence of foundation spacing and soil behaviour from the confounding effects of absolute strength, allowing for a more precise assessment of how spatial interactions modify the expected performance.

Figure 5

Normalised analytical and numerical results for deterministic values of soil parameters for each approach: (a) ideal cohesive, (b) purely frictional, and (c) TBC.

As a result of normalising the outcomes, an increase in BC of up to 15% is observed for foundations spaced 0.5 m apart in Approach 2. In Approach 3, the increase is approximately 10%, which is 5% higher than the BC of a single foundation with a weight of 2 m, whereas in Approach 1, it remains close to 1%. This enhancement does not evolve uniformly with increasing spacing. In Approach 3, the strengthening effect diminishes almost linearly with increasing distance. In contrast, in Approach 2, the BC decreases until the spacing reaches approximately 2.5B, after which it rises again to 8% for larger foundation spacing. This effect is not readily apparent and warrants further analysis. This comparative perspective underscores that the magnitude and even the direction of interaction effects are not universal but depend strongly on the underlying soil mechanics. For engineering design, this means that the potential benefits or drawbacks of placing foundations in proximity cannot be inferred from one soil type to another without careful consideration of the specific context.

The results of probabilistic analyses performed using RFEM are shown in Figures 6 and 7. Figure 6 illustrates the mean BC for probabilistically modelled soil parameters across three Approaches: (1) perfectly cohesive, (2) purely frictional, and (3) TBC. The results reveal distinct patterns for each case. The results in Figure 6 are normalised under the same conditions as in the previous subsection, by dividing the adequate analytical value from Table 1. In Approach 1 (Figure 6a), the mean BC is approximately 20% lower than the analytical value, indicating a significant reduction due to spatial variability in cohesion. For Approach 3 (Figure 6c), the mean values are also lower, but the reduction does not exceed 10% of the analytical mean for B = 1 m. Notably, when the gap between foundations is 1 m or less, the mean BC increases above the nominal value.

Figure 6

Mean BC for probabilistic values of soil parameters for each approach: (a) ideal cohesive, (b) purely frictional, and (c) TBC.

Figure 7

BC COV for probabilistic values of soil parameters for each approach: (a) perfectly cohesive, (b) purely frictional, and (c) TBC.

The non-cohesive (purely frictional) soil (Figure 6b) exhibits the most pronounced effect of spatial variability: for small gaps (0.25–1.0 m), the mean BC almost mirrors the deterministic results, while for gaps greater than 1 m, the mean values are slightly lower than their deterministic counterparts. The influence of the SOF _x on the mean BC is minimal, as evidenced in Figure 6.

Figure 7 shows the coefficient of variation (COV) for each Approach, which indicates a stronger relationship with SOFs. For SOF _x = 1 m, the COV is consistently lower, regardless of the soil type. At larger scales (SOF _x = 10 and 100 m), the COV values are similar, with variability generally increasing with scale. Furthermore, the COV exhibits fluctuations for gaps less than 1 m, independent of the soil model adopted.

A further point of interest is the dependence of COV on soil type. The highest COV is observed in the perfectly cohesive soil (Figure 7a), while the purely frictional soil (Figure 7b) displays the smallest COV. The natural soil (TBC, Figure 7c) presents intermediate values. This pattern suggests that cohesion variability directly influences the variability in the BC of adjacent foundations.

To synthesise the effect of soil type on the observed results, the author introduces a global Stuart reinforcement index ( ζ G ) {\zeta }_{{\rm{G}}}) , defined as the ratio of the mean BC from both foundations to the analytical BC for a single footing (equation 6). (6) ζ G = Mean BC from both foundations Analytical BC value for a single footing {\zeta }_{{\rm{G}}}=\frac{\text{Mean}\text{BC}\text{from both foundations}}{\text{Analytical BC value for a single footing}}

Figure 8 presents the global efficiency factor values for all considered approaches. The analytical factor suggested by Stuart (1963) was in the range of 1–2.6 with axial spacing Δ from 4B to 5B (GAP 5B–6B) for purely frictional soils (as Approach 2), while for cohesive soils (like Approach 1), the range of factor was 1–1.35 for Δ = 3B (GAP = 4B). In the current analysis, bearing-capacity reinforcement is evident for GAPs less than 2B; in the probabilistic context, the magnitude of this reinforcement for natural soil under Approach 3 does not exceed 10%.

Figure 8

Performance efficiency factor for all approaches (A1 – Approach 1 – perfectly cohesive; A2 – Approach 2 – purely frictional; A3 – Approach 3 – TBC).

In Approach 2, the gain is most pronounced, reaching up to 20%. Conversely, cohesive soil demonstrates a reduction in the coefficient of up to 20%, regardless of the mutual position of the foundations.

The probability of failure for adjacent shallow foundations, as illustrated in Figure 9, is evaluated using equation 5 under two principal design assumptions regarding the threshold values. First, with the values of q a , d {q}_{a,d} , and second, q a , d {q}_{a,d}\hspace{0.25em} multiplied by the Stuart coefficient (ζ, Figure 8). Each row in Figure 9 corresponds to a specific soil type: Approach 1 (Figure 9a and b), Approach 2 (Figure 9c and d), and Approach 3 (Figure 9e and f). Each subfigure in Figure 9 depicts the relationship between probability of failure and GAP for two sets of results: values based on the analytical BC divided only by FS, represented by lines with markers, and values where the analytical BC is divided by FS and further multiplied by the Stuart coefficient obtained for each approach, represented by markers only. In each set, three values of SOF were considered.

Figure 9

Probability of failure. (a and b) Approach 1 – perfectly cohesive soil, (c and d) Approach 2 – purely frictional soil, and (e and f) Approach 3 – cohesive-frictional soil (TBC).

A previous study by the author [25] demonstrated that the probability of failure strongly depends on the selected FS value, and that this safety level directly affects the estimation of the reliability index ( β ) \left(\beta ) . As illustrated in the Figure 9a, c and e the probability of failure exceeds 5% when FS = 1. That situation corresponds to a reliability index beta of around 0.5. To achieve a failure probability of p f ∼ 10 − 4 {p}_{{\rm{f}}}\sim {10}^{-4} , as recommended in Eurocode 7, the reliability index must exceed β ≥ 3.55 \beta \ge 3.55 . This requirement necessitates adopting different FS values for each soil type and approach, reflecting their inherent variability and associated reliability demands, as represented by column (2) in Figure 9b, d, and f.

Under the baseline scenario (FS = 1), a strong correlation is observed between failure probability and mean BC (as shown in Figure 6). When the Stuart coefficient is applied, the probability of failure stabilises at approximately 35% for all soil types, regardless of the correlation structure or relative foundation position.

For Approach 1 – perfectly cohesive soils (Figure 9b), achieving the target reliability index β ≥ 3.55 required a FS = 4.5. In this case, the probability of failure approaches zero when the Stuart coefficient is included, indicating a highly conservative design. However, if the limit value is based solely on the standard analytical capacity (without the Stuart adjustment), the influence of spatial correlation becomes more pronounced, underscoring the importance of accounting for soil variability in reliability assessments.

In contrast, for Approach 2 (Figure 9d), the required SF to achieve β ≥ 3.55 is significantly lower (FS = 2.2). Here the effect of the Stuart coefficient is less pronounced, whereas the impact of the correlation structure is more pronounced. As the scale of fluctuation increases, so does the probability of failure in both scenarios. The highest failure probabilities are observed for GAPs between 1 and 3 m, indicating that interaction effects are most significant within this range under variable soil conditions.

For Approach 3 (Figure 9f), the effect of the correlation structure dominates. With a small fluctuation scale (SOF = 1 m), the probability of failure remains close to zero in both scenarios. As SOF increases, failure probability rises, particularly when the Stuart coefficient is applied, highlighting the soil’s sensitivity to spatial variability and interaction effects. A particularly noteworthy observation is the influence of GAP on failure probability. Across all soil types, the failure probability increases sharply for GAP values between 0 and 3 m, reflecting strong interaction effects at close distances. For spacings greater than 3 m, the probability of failure stabilises or decreases slightly. This indicates that foundations begin to behave independently, and the interaction effect diminishes. Such a trend reinforces earlier findings that adjacent footings exhibit significant load-bearing interaction when placed in proximity; however, beyond a particular spacing, this effect wanes, and the probability of failure approaches that of isolated foundations.