3D bearing capacity probabilistic analyses of footings on spatially variable c–φ soil

Abstract

The paper examines three-dimensional (3D) analyses of the load bearing capacity of square and strip footings on a spatially variable cohesive–frictional (c–φ) soil. Embedded, rigid and smooth footings are considered. Cohesion and friction angle are modelled using random fields. The vertical scale of fluctuation (SOF) is based on earlier studies of natural soils. For the square footing problem, the horizontal SOF is the subject of a parametric study. For the strip footing, only one value of the horizontal SOF is assumed, and the parametric study concerns the modelled length of the foundation. To generate the fields of strength parameters, the Fourier series method is employed. For both problems analyzed, the influence of the parameters on the probabilistic characteristics of the bearing capacity of the soil is investigated. For the square footing problem, both the mean value of bearing capacity and its coefficient of variation increase with an increasing horizontal fluctuation scale, but they stabilize at higher values of SOF. For the strip footing problem, it appears that 2D modelling is a very conservative approach, and for the precise assessment of the probability of foundation failure, the footing should be modelled in 3D.

1 Introduction

The past 20 years have seen increasingly rapid advances in the field of random analysis of foundations based on modelling spatial variability in soil properties using random fields. Many researchers have studied this problem [1, 4, 9, 14, 16, 20, 29, 32, 36]. The pioneering work on this subject by Griffiths and Fenton [15] and Fenton and Griffiths [12] analyzed the random bearing capacity of strip foundations. The academic problems were initially limited to embedded foundations on purely cohesive soils modelled with isotropic fields, but have recently been developed for studies that involve analyses of design situations involving foundations on frictional as well as cohesive soils. For example, in the work of Pieczyńska-Kozłowska et al. [30], the bearing capacity of an embedded strip foundation on a natural cohesive–frictional soil was analyzed.

The random analysis of such nonlinear problems requires significant computing power. This is probably the reason why papers examining the three-dimensional (3D) effects of geotechnical problems with variable soil are limited in number. One exception was the early work of Fenton and Griffiths [13] in which foundation settlement was analyzed in 3D. However, the problem in that case was linear, which strongly simplified the computation process. The plasticity problem for spatially variable soil in 3D started to be analyzed during the previous decade, and until recently [especially for a significant number of Monte Carlo simulations (MCS)], was beyond the computing capabilities of personal computers. A good example of such an analysis was the series of papers presented by Hicks et al. [17,18,19, 26, 27]. To analyze a sufficient number of realizations of the large-scale 3D slopes, the authors employed high-performance computing clusters.

The first probabilistic 3D analysis of a square footing on soil with strength modelled with random fields was carried out by Al-Bittar and Soubra [3] (and in the Ph.D. thesis of Al-Bittar [2]). The FLAC3D finite difference code was used to solve the boundary value problem. Because the direct use of the finite difference code was computationally demanding, a metamodel was developed based on sparse polynomial chaos expansion combined with global sensitivity analysis. Simoes et al. [35] used an upper bound formulation of finite element limit analysis (FELA) in order to check the influence of modelling strip foundations on spatially variable 3D soil with regard to the probability characteristics of structure bearing capacity. For these simulations, the Latin hypercube sampling technique was used, which allowed the limiting of a satisfactory number of simulations compared to the crude Monte Carlo technique. The computations were performed using a cluster with a total number of 102 cores. Both of the above studies concerned undrained (or purely cohesive) soil media and not embedded foundations.

The surface square footing was also analyzed in one of the earlier studies by the authors of the present paper [23]. The local average subdivision (LAS) procedure was utilized as the method for random field generation. The analysis involved both purely cohesive and cohesive–frictional soils. The boundary value problem was solved using the FLAC3D code. Since this work was a preliminary study, it was simplified by assuming that soil variability takes place in only the vertical direction. (The soil properties were assumed not to vary in the horizontal direction so that the horizontal scale of fluctuation was infinity.)

It should be noted that the task of analyzing the random bearing capacity of foundations in 3D is not only computationally challenging (as evidenced by few publications on the topic), but also very important. Rectangular footings in industrial structures (i.e. where heavy loads can occur and the risk of failure is worth assessing) is the most common foundation type. The bearing capacity of a rectangular footing is, of itself, a 3D problem. The 3D modelling of a strip foundation located in spatially variable soil is also justified. As shown in [35], the modelling of a strip foundation in 2D can lead to much higher values of bearing capacity coefficients of variation than would result from a full 3D analysis. Thus, a 2D analysis of a strip footing may be a too conservative approach, which is of great importance when evaluating the probability of foundation failure.

The present work, which is an extension of earlier studies, aims to analyze the random bearing capacity of embedded square and strip footings resting on a natural spatially variable cohesive–frictional soil. FLAC3D is used to solve both types of boundary value problem. The well-recognized soil known as Taranto clay is chosen as a cohesive–frictional soil. Since no information about the horizontal scale of fluctuation for that soil can be found in the literature, a parametric study of that value is performed to address the square footing problem. Random fields are generated using the Fourier series method (FSM) and then transformed to match the distributions of cohesion and friction angle. The same methodology is used to analyze the influence of 3D modelling of strip footings on the probabilistic bearing capacity characteristics.

This study is structured as follows. The next section presents the basic assumptions regarding the probabilistic modelling method, adopted distributions and values of the modelled soil parameters. The numerical modelling is presented in the subsequent section followed by the presentation of the obtained probabilistic results. The work is summarized by a section discussing the conclusions of the study.

2 Basic assumptions

2.1 Random field modelling

Modelling the spatial variability in soil properties by random fields seems to be a reasonable solution if a random bearing capacity is under consideration, and it is well described in the literature [12, 23, 30, 39]. In the present study, two soil properties are modelled using a random field: the cohesion and the angle of internal friction. The random fields under consideration are assumed to be weakly stationary and ergodic [14]. For both of them, the same correlation structure has been assumed by applying the following correlation function (often called the Markov-type function):

$$\rho \left( {\tau_{x} ,\tau_{y} ,\tau_{z} } \right) = { \exp }\left( {\frac{{ - 2\left| {\tau_{x} } \right|}}{{\theta_{x} }} + \frac{{ - 2\left| {\tau_{y} } \right|}}{{\theta_{y} }} + \frac{{ - 2\left| {\tau_{z} } \right|}}{{\theta_{z} }}} \right)$$

(1)

where (τ_x, τ_y, τ_z) = (x₂-x₁, y₂–y₁, z₂–z₁). Parameters θ_x, θ_y and θ_z are scales of fluctuation (correlation radii) in the directions x, y and z, respectively. In the following, the acronym SOF is used instead of the term “scale of fluctuation.”

Function (1) is commonly used in soil properties modelling. Accepting the same correlation function (with identical values for the SOF) for both the friction angle and cohesion random fields suggests that the spatial correlation of properties in a soil is an effect of its geological history rather than a “specific character” of a given soil parameter. The problem of SOF estimation has been intensively explored in the literature [5, 6, 21, 31, 33, 38, 40]. Generally, methods are well established for the vertical SOF (“variability with depth”), but not satisfactorily developed for the horizontal SOF. For this reason, the effect of horizontal SOF on bearing capacity probabilistic characteristics in this study will be the subject of a parametric study. In earlier papers, evaluating bearing capacity [12, 15], a simplification of the correlation structures has been imposed by assuming that the SOF is identical in all directions. On the other hand, as indicated by Cherubini [10], the value of the horizontal scale of fluctuation is usually greater than the vertical one. It was demonstrated in papers [30, 39] that, for a given vertical SOF, using the same value of SOF in both the vertical and horizontal directions may lead to smaller values of bearing capacity coefficient of variation compared to cases where the horizontal SOF was significantly greater than the vertical one. For this reason, the present study assumes significantly different values of SOF in the vertical and horizontal directions.

The statistical sampling carried out in the present study utilizes the Monte Carlo technique. Therefore, the numerical simulation of a random field under consideration plays a vital role. Here, based on one of the pioneering studies by Fenton and Griffiths [12], random fields that model soil strength parameters are obtained by the appropriate transformation of normal random fields (i.e. random fields with normal finite-dimensional distributions). For example (as suggested in papers by Griffiths and Fenton [15] and Fenton and Griffiths [12]), if the cohesion is assumed to be a log-normal random field Y(x) (where x ∈ D, where D is a three-dimensional domain), the transformed field is obtained from a normal random field X(x) by applying the following transformation:

$${Y}\left( {\mathbf{x}} \right) = {\exp}\left( {{X}\left( {\mathbf{x}} \right)} \right)$$

(2)

At a specified point x, the probability density function (p.d.f.) of Y(x) takes the form

$$f\left( {t} \right) = \left\{ {\begin{array}{ll} {\frac{1}{{{t}\sigma _{{\ln {Y}}} \sqrt {2{\pi}} }}\exp \left\{ { - \frac{1}{2}\left( {\frac{{\ln {t} - {\mu} _{{\ln {Y}}} }}{{{\sigma} _{{\ln {Y}}} }}} \right)^{2} } \right\},} \hfill & {{t} > 0} \hfill \\ 0, \hfill & {{t} \le 0} \hfill \\ \end{array} } \right\}$$

(3)

where parameters μ_lnY and σ_lnY are the mean value and the standard deviation of the underlying normal distribution X(x), respectively. The mean value μ_Y and the standard deviation σ_Y of the log-normal distribution Y(x) are related to μ_lnY and σ_lnY by

$$E\left[ Y \right] = \mu_{Y} = \exp \left( {\mu_{\ln Y} + \frac{1}{2}\sigma_{\ln Y}^{2} } \right)$$

(4)

$$\sigma_{Y}^{2} = \mu_{Y}^{2} \left( {{ \exp }\left( {\sigma_{{{ \ln }Y}}^{2} } \right) - 1} \right)$$

(5)

With regard to the modelling of the friction angle, due to the bounded range of this parameter, the probability distribution function should be a function of a bounded support, where the support of a function is defined by: supp(f) = cl{x: f(x) ≠ 0} (cl denotes the closure operation). In order to meet this requirement, Fenton and Griffiths [12, 14] suggested the following transformation:

$$Y\left( {\mathbf{x}} \right) = a + \frac{1}{2}\left( {b - a} \right)\left[ {1 + { \tanh }\left( {\frac{{m + sX_{0} \left( {\mathbf{x}} \right)}}{2\pi }} \right)} \right]$$

(6)

where X₀(x) is the standard normal random variable. The support of the p.d.f. of Y(x) is the interval [a,b], and m and s are the distribution’s parameters. The p.d.f. of Y(x) takes the form:

$$f\left( t \right) = \left\{ {\begin{array}{*{20}l} {\frac{{\sqrt \pi \left( {b - a} \right)}}{{\sqrt 2 {\text{s}}\left( {t - a} \right)\left( {b - t} \right)}}{ \exp }\left\{ { - \frac{1}{{2{\text{s}}^{2} }}\left[ {\pi { \ln }\left( {\frac{t - a}{b - t}} \right) - m} \right]^{2} } \right\},} \hfill & {t \in \left( {a,b} \right)} \hfill \\ 0 \hfill & {t \notin \left( {a,b} \right)} \hfill \\ \end{array} } \right.$$

(7)

If m = 0, the distribution is symmetric with respect to the centre of the variability interval [a,b]. The parameter s is associated with the standard deviation of Y(x); however, there is no closed form relation between them. By means of the third-order Taylor’s expansion, the following approximation can be obtained:

$$\sigma_{Y}^{2} \approx \frac{1}{4}\left( {b - a} \right)^{2} \frac{{s^{2} }}{{4\pi^{2} + s^{2} }}$$

(8)

2.2 The soil characterization: Taranto clay

As mentioned in the introduction, the actual problem considered in this paper is the analysis of square and strip footings on a cohesive–frictional soil. To account for the physical properties of real soil, the present analysis is based on extensive research (including statistical measures and the determination of a vertical scale of fluctuation) of natural deposits of Taranto clays, which was carried out by Cafaro et al. [7, 8]. These deposits are situated in southeast Italy (Fig. 1). The Taranto clays are stiff and over-consolidated and are mainly composed of illite and kaolinite minerals.

Fig. 1

Geographical location of Taranto clays

The statistical data in the cited studies were collected by means of CPT testing in 15 boreholes. The testing demonstrated the presence of two clay layers: an upper layer composed of yellow clay and a lower layer composed of blue–grey clay (Fig. 2). Both clays can be classified as clayey silt to silty clay.

Fig. 2

Soil layers recognized by CPT tests

The variability in the strength parameters for the Taranto clays was determined by means of the statistical treatment of the CPT results. Data from five boreholes (selected from the fifteen boreholes presented in Fig. 2) for both the yellow and the blue–grey clays were analyzed. In order to obtain the residuals, a low-order polynomial trend (no higher than a quadratic) was removed by regressing the cone resistance (q_c) values using the ordinary least squares method. Furthermore, the geostatistical boundaries of the two Taranto clays were investigated by means of the interclass correlation coefficient [40]. The results showed the presence of a lithological boundary at a depth of approximately 11 m. The mean trend in q_c with depth z, given by the trend equation f(z) as well as the vertical scales of fluctuation that were obtained, are presented in Table 1. Unfortunately, the horizontal fluctuation scale was not investigated.

Table 1 Trends and vertical scales of the fluctuation from CPTs within upper and lower clay

The scale of fluctuation was evaluated after the removal of the mean trend line from the q_c readings. As regards strength parameters, two approaches are possible. The first one assumes evaluation of the strength parameters (φ and c) based on q_c readings, followed by estimation of the mean values μ_φ and μ_c as well as standard deviations σ_φ and σ_c. The values determined for μ_φ, σ_φ and μ_c, σ_c are used as the first two statistical moments of the random field of φ and the random field of c, respectively. The above approach allows the weak stationarity of both random fields under consideration to be preserved. The second approach allows the use of trends (of mean value and standard deviation) of φ and c, by adding them during the simulation of random fields. However, in such an approach, the fields of φ and c are no longer stationary. Moreover, as was demonstrated in [24], the friction angle displayed almost no trend despite the evident trend in q_c. In the present study, the first approach was utilized.

As described in Sect. 2.1, suggested distributions for the cohesion and friction angle are log-normal [Eq. (3)] and bounded [Eq. (7)] distributions, respectively. Based on the collected data and by applying Eqs. (4), (5) and (8), the parameters of the above-mentioned distributions were identified. The parameters of the probability distributions of soil properties obtained from the Taranto blue–grey clay measurements that were utilized in the present study are presented in Table 2.

Table 2 Random mechanical properties characterization of blue–grey clay

The same parameters as presented in Table 2 were applied in study [39] following the suggestion of Cherubini [11]. From the data given in Table 2, the relatively high value of coefficient of variation (CoV) of about 0.55 is obtained for c. This could be a consequence of applying the first approach mentioned above. However, even higher values of the coefficient of variation in c have been reported in some earlier papers, e.g. [25]. Also, some investigation of the clay deposits of Lower Silesia (the region where the authors of the present study live) shows high values of the COVs of c as reported in the Ph.D. thesis of Thao [37]. Additionally, using the same values of parameters as in earlier work of the authors (e.g. [30]) allows the effects of 3D and 2D modelling to be compared for similar variability in soil parameters.

Due to extensive computational time, only one value of the vertical fluctuation scale, θ_z = 0.5 m, was selected as a representative value for further computations.

3 Numerical model

The 3D probabilistic analysis presented in this paper concerns the bearing capacity of square and strip footings. The numerical computations were performed in the FLAC3D environment. The elastic-perfectly plastic Mohr–Coulomb material (with non-associated flow rule) was employed as a soil model. The Young’s modulus, Poisson’s ratio, unit weight and dilation angle were given non-random values of E = 35 MPa, ν = 0.35, γ = 19 kN/m³ and ψ = 0°, respectively. Only the parameters for the strength criterion (namely cohesion c and friction angle φ) were assumed to be random. These parameters were modelled using random fields to recreate the probabilistic characteristics of Taranto clays (i.e. mean values, standard deviations, distribution types and vertical SOF as presented in Table 2). The fields of c and φ were assumed to be stochastically independent.

As mentioned in Sect. 2.1, the exponential function (1) has been chosen to establish the correlation structure in random fields within this study. The function (1) is positively defined (a sufficient and necessary condition according to Bochner’s theorem [34]), reflects separability of the fields under consideration, gives the possibility of distinct scales of fluctuation in three directions and is suggested by many researchers (e.g. [28]).

To generate random fields with the autocorrelation function of type (1) for a given mesh, the algorithm proposed by Jha and Ching [22] (FSM) based on the Fourier series expansion was employed. The Fourier series method (like most of the generation methods) generates normally distributed fields. Because of that, for the log-normal field of cohesion, the underlying normal distribution parameters were first found by solving Eqs. (4) and (5) for μ_lnY and σ_lnY. Then, the generation was carried out, and the transformation (2) then applied to all generated values. In the case of the internal friction angle field, the parameter s = 2.27 [see Eq. (8)] was evaluated to obtain a measured standard deviation for the given range of bounds. Next, the generated standard normal field was transformed using Eq. (6). For both problems analyzed (square footing and strip footing), the vertical SOF was assumed to be 0.5 m. Since the value of the horizontal scale of fluctuation for Taranto clay was not investigated, some typical values of horizontal SOF have been assumed as specified in the next subsections. In all the present analyses, both horizontal SOFs (namely θ_x and θ_y) were assumed to be equal. Example realizations of the log-normal field of cohesion and the field of friction angle [with p.d.f. defined by Eq. (7)] for the mesh applied to the square footing problem are presented in Fig. 3. Both the range of values for the typical realization as well as the size of strongly correlated clusters for the given SOFs are shown.

Fig. 3

Contour diagrams of example realizations of random parameter fields generated using FSM for 8 m × 8 m × 3.5 m domain, θ_x = θ_y = 1 m and θ_z = 0.5 m: a the field of cohesion (log-normal) and b the field of the friction angle [distribution with p.d.f. given by Eq. (7)]

3.1 Square footing problem

For the square footing problem, the footing was assumed to be rigid, smooth and only vertically loaded. The horizontal dimensions of the footing were assumed to be 1.17 × 1.17 m. The boundary conditions were set as follows: The displacement of nodes lying on the external x–z and y–z planes (where z is the vertical direction) of the cuboid domain was fixed in directions normal to the planes and on the bottom x–y plane in all three directions. The 1-m embedment was simulated using the additional loading of 19 kPa applied to the upper surface of the domain everywhere except the foundation area. The footing was modelled as 7 × 7 nodes spaced 0.17 m (in both directions) with an applied constant velocity of displacement of the order 5 × 10⁻⁵ m/time step. The bearing capacity was defined as the maximum vertical force which can be applied to the footing. Knowing this value, the value of the maximum (average) vertical stress under the footing can also be easily calculated (by dividing the force by the area of the footing). The force was calculated as the sum of vertical reactions in all the footing’s nodes. This value stabilized after several thousand steps when it reached the limit value of the load bearing capacity.

The dimension of the domain (surrounding the footing) was found to be large enough to minimize the effects of the boundary on the obtained bearing capacity. To determine that dimension, a non-random study was first conducted. The c and φ values were assumed to be equal to their mean values, i.e. 36 kPa and 20°, respectively. All other parameters were assumed to be the same as specified in the previous section. Due to the symmetry of the non-random problem, only one quadrant of the problem domain was considered. The obtained values for the footing bearing capacity Q_f for the dimensions of the domain’s quadrant equalling 3, 2 and 1 m, together with respective images of the plastic zones, are presented in Fig. 4. The blue, brown and pink colours in the figure denote that the cell is in an elastic state (without plastic history), plastic state (yielding in shearing) or elastic state with plastic history, respectively. The nodes associated with the footing are marked with vertical black lines. As shown in Fig. 4, the boundary effect is limited. The case “a” in which the whole plastic zone is contained in the domain is treated as the reference case. It was verified that any enlargement in the domain dimension did not change the obtained value for the bearing capacity. Compared to situations in which the plastic zone only “touches” the domain’s boundary (case “b”), no significant decrease in load bearing capacity was observed. Only when the area of contact between plastic zones and the boundary was very large (case “c”), there was a significant decrease in the bearing capacity observed. Thus, if all or almost all of the plastic zones were contained in the domain, the influence of the boundary was negligible.

Fig. 4

Plastic zones under the square footing. The dimensions of the analyzed part (domain quadrant) was assumed to be a 3 m, b 2 m or c 1 m

A similar analysis was conducted for the problem with strength parameters modelled with random fields. Obviously, in that case, the problem was no longer symmetrical; the whole domain needed to be analyzed instead of just one quadrant. For a few assumed horizontal SOFs ranging from 1 to 30 m, the influence of the domain size on the mean value of bearing capacity was tested. Since it seemed that the cases of larger SOF were more sensitive to the domain being too small, this influence was carefully investigated in the case of a horizontal SOF of 30 m. Finally, the domain size of 8 × 8 × 3.5 m was found to be satisfactory, and enlarging the domain did not change the mean value by more than 0.5% (although that difference could be the effect of probabilistic modelling). With the assumed uniform grid with finite difference zone dimensions of 0.17 × 0.17 × 0.17 m, a mesh of 48 × 48 × 21 zones was obtained. In Fig. 5, some exemplary plastic zones obtained in different realizations for that domain are presented. The same colour coding as shown in Fig. 4 is used. Although the entire domain was analyzed, only half of it is shown in the figure. This is to show the plastic zones are not only on the boundary of the domain, but also directly under the footing. Two of the presented realizations were obtained for minimal (1 m) and the other two for maximal (30 m) considered horizontal SOFs. (The zones for SOF = 30 m are easy to recognize as they are much more symmetrical.) As shown in Fig. 5, in each case, almost the whole plastic zone is contained inside the domain. However, it should be noted that for the horizontal SOF of 30 m, realizations in which the plastic zone is “touching” the boundary are much more frequent. It was also found that for the cases of higher horizontal SOF, reducing the domain dimension resulted in numerical errors that, even though they were few in number, could significantly influence the mean value of the bearing capacity.

Fig. 5

Representative images of plastic zones for half of the considered domain: a, b for θ_x = θ_y = 1 m; c, d for θ_x = θ_y = 30 m

For the given domain, the influence of the horizontal SOF value on the probabilistic characteristics of bearing capacity was investigated. Five different horizontal SOFs (i.e. 1, 2, 5, 10 and 30 m) have been tested. For each case, 300 Monte Carlo simulations (MCS) have been computed. This number was limited due to the calculation time: obtaining 300 solutions took about 2 weeks of work using a single processor core on a modern PC.

3.2 Strip footing problem

The same procedure for modelling was also used to analyze the impact of the 3D field effect on the bearing capacity of the strip foundation. The 1.17 m wide perfectly smooth and rigid strip footing embedded 1 m into cohesive–frictional soil was considered. All the soil parameters as well as the field generation process and the procedure of finding the critical load were the same as described previously. This time, the bearing capacity was defined as the maximum allowable vertical force per metre of the foundation length. The horizontal dimension of the domain in the direction perpendicular to the strip footing and the vertical dimension were once again chosen to avoid boundary effects. Based on the analysis of the influence of the domain dimensions on the mean value of the bearing capacity, a width (perpendicular horizontal dimension) of 14 m and depth (vertical dimension) of 6 m were found to be satisfactory. For these domain dimensions, plastic zones quite often “touched” the boundary, but simulations carried out with increased dimensions hardly affected the results (apart from a significant extension of the calculation time). In fact, the obtained difference in mean value was of the order of 0.2%. Because the uniform grid of cubic zones with all dimensions equal to 0.17 m was used (identically to the square footing problem), the number of zones in a single slice of the grid (in the plane perpendicular to the direction of the strip footing) was 84 × 36.

In the case of a strip footing, the impact of modelling the foundation in 3D (on a 3D spatially variable soil) on the footing load carrying capacity results was investigated. A parametric study of that case did not concern the horizontal SOF, but investigated the effect of the modelled length of the strip footing. Only one case of horizontal SOF was analyzed, namely θ_x = θ_y = 2 m. The analysis included five series of calculations, each for different lengths of strip footing, i.e. 0.17, 1, 2, 4 and 6 m (1 zone, 6 zones, 12 zones, 24 zones and 36 zones, respectively). The first series (0.17 m) was considered to be a 2D problem, but in the entire series, there was a general assumption that only part of the foundation was being modelled because the foundation runs through the entire modelled length of the domain, and neither of its ends were modelled (see Fig. 6).

Fig. 6

Representative plastic zones obtained for four different lengths of modelled strip footing (0.17 m, 1 m, 2 m and 4 m). The footing nodes are denoted with vertical black lines

Similarly, as in the case of the square footing, each series consisted of 300 MCS. In this case, the analysis was even more time consuming. In particular, 300 MCSs of a 6-m-long portion of the strip footing required about 6 weeks of computational work for a single core. Representative plastic zones (using the same colour coding as previously) obtained for four different modelled lengths of footing are presented in Fig. 6. (The results of the fifth series, which are not presented here, were qualitatively similar to series 4.) The results obtained for the bearing capacity of both the square and strip footings are presented in the next section.

4 Results and discussion of the probabilistic analysis

4.1 Bearing capacity of square footing

As mentioned in the previous section, the numerical simulation for the square footing problem was conducted to determine the influence of the horizontal SOF on the probabilistic results. For each considered SOF, N = 300 realizations of MCS were carried out. Based on the collected data, in each case, both the mean and coefficient of variation (CoV) of Q_f were estimated. The diagrams showing the influence of these two characteristics as a function of the horizontal SOF are presented in Fig. 7. For comparative purposes also, the mean value of Q_f obtained for the non-random case (where the values of c and φ were assumed to be their mean values, i.e. 36 kPa and 20°, respectively) is presented in Fig. 7a. As shown in the figure, for increasing values of the horizontal SOF, both the mean value and the CoV of bearing capacity increase, but these soon stabilize for higher values of SOF. Slight deviations from the observed trend obtained in both diagrams for horizontal SOF between 10 and 30 m may have resulted from the limited number of realizations. Since all mean values obtained lie in a very small range (about 2–3% difference between all obtained values), the observed increase in the mean is very small. On the other hand, the change in CoV is quite significant. Its value increases from 0.15 (obtained for the horizontal scale of 1 m) to 0.36 (obtained for the horizontal scales of 10 m or 30 m). It is worth noting that the maximum value of CoV obtained for high values of SOF as well as the general trend in the CoV function is in good agreement with the results obtained in the previous study of a strip foundation modelled in 2D on the same soil for high values of horizontal SOF [30].

Fig. 7

Influence of the horizontal fluctuation scale on a mean value and b CoV of square footing bearing capacity

The probability distributions of the bearing capacity for the individual series of calculations were estimated and tested using Kolmogorov’s goodness-of-fit test. For the hypothetical log-normal distribution of bearing capacity, the test results were in the range of 30–80%. The obtained estimated distributions for different values of the horizontal scale of fluctuation are presented in Fig. 8.

Fig. 8

Estimated distributions of the bearing capacity Q_f for different values of horizontal SOF

Since the sample size was quite limited (N = 300), the estimated distributions presented in Fig. 8 should be treated as only a rough approximation, especially in cases when the tails of the distributions have not been precisely estimated. Nevertheless, a simple reliability analysis was performed by assuming that the obtained distributions were the correct ones.

Let us define the probability of failure as the probability that the bearing capacity Q_f is lower than a certain allowable (minimal) value Q_d, namely

$$p_{\text{f}} = p\left( {Q_{\text{f}} \le Q_{\text{d}} } \right) = \varPhi \left( { - \beta } \right),$$

(9)

where Φ is the cumulative distribution function of the one-dimensional standard normal distribution and β is the reliability index. According to the Eurocode regulations, the target reliability index for the ultimate limit states of RC2 reliability class structure (ordinary buildings and their foundations) having a design working life of 50 years has been specified as 3.8. Using this value, the values of Q_d as a function of horizontal SOF were calculated.

The obtained results are presented in Fig. 9. As shown, in the case of small values of SOF, the allowable value of Q_d is significantly higher. In contrast, for higher values of SOF, which seem to be the situation for the natural soil [10], the value is much lower but quite stable. Figure 9 corresponds well to Fig. 7. Both the mean value and standard deviation are almost constant (see Fig. 7) for SOF greater than 10, resulting in an almost constant probability of exceeding the allowable value Q_d.

Fig. 9

Allowable value Q_d of square footing bearing capacity corresponding to β = 3.8 versus the horizontal SOF

4.2 Bearing capacity of strip footing

The analysis of the bearing capacity of the strip footing aimed to determine the impact of 3D modelling of the footing and soil on the obtained probabilistic results. As previously mentioned, the parametric study concerned the modelled length of the foundation. Only one horizontal SOF value (i.e. θ_x = θ_y = 2 m) was assumed. For each of the foundation lengths considered, N = 300 MCSs were performed.

The estimated means and CoVs are presented in Fig. 10. As previously mentioned, for comparative purposes the value of Q_f obtained for the non-random case (where the values of c and φ were assumed to be their mean values, i.e. 36 kPa and 20°, respectively) is also presented in the figure. The minimum length of the modelled foundation equals 0.17 m, corresponding to a single slice of the finite difference mesh as indicated in the previous section. The increase in the mean value of the bearing capacity, which occurred for higher values of the modelled length of the strip foundation, was not significant; the maximum increment of the mean value was 6%. In contrast, the CoV of the bearing capacity significantly decreased for higher values of the model length. The changes in both the mean and the CoV of bearing capacity caused by the increase in the model length were the effect of averaging the bearing capacity along the foundation.

Fig. 10

Influence of the modelled length of the foundation on the probability characteristics of its bearing capacity Q_f : a mean value and b CoV for horizontal SOF of 2 m and vertical SOF of 0.5 m

Since the range of values of mean bearing capacity is limited and the graph becomes flat for larger modelled lengths of the foundation, it seems that the value obtained for the length of 6 m is close to an “asymptotic” value. On the other hand, the diagram of the CoV is much steeper which indicates that for very long lengths of the foundation, the CoV of the bearing capacity can approach a much smaller value than the 0.1 obtained for 6 m. These results seem to agree with the previously mentioned averaging effect. It should be noted that these results were obtained for a horizontal SOF equal to 2 m, which is a relatively low value. Up until now (due to high computation costs), we have no results for higher values of the horizontal SOF in the case of the strip foundation. However, one can expect that for the horizontal SOF equal to infinity, there should be no influence of the modelled length of the foundation on the results.

Once again, in spite of the limited sample size, the probability distributions of the bearing capacity of the strip foundation were estimated and tested using the goodness-of-fit test. The results obtained for the Kolmogorov test were between 30 and 70%, and the hypotheses of log-normal distributions were not rejected. The hypothetical distributions for differently assumed lengths of the model are presented in Fig. 11.

Fig. 11

Hypothetical distributions of bearing capacity obtained for different lengths of strip foundation

Once again, a simple reliability assessment was carried out. The allowable values of the bearing capacity Q_d corresponding to the value of the reliability index β = 3.8 were evaluated. The results are presented in Fig. 12.

Fig. 12

Allowable values Q_d of strip footing bearing capacity corresponding to β = 3.8 versus the modelled length of footing l

As shown in Fig. 12, the allowable value Q_d increases with the longer modelled length of a foundation. The change is significant: The allowable value Q_d for the modelled length of 6 m at 470 kN/m is more than 1.5 times greater than the Q_d of 270 kN/m obtained for 0.17 m. It should be noted that the case of the minimum model length (equal to the length of the single zone) corresponds to the 2D random analysis of the strip foundation. The obtained results indicate that the 2D reliability analysis can be significantly conservative in comparison with 3D analysis. It seems that modelling such structures in 3D is necessary for the precise assessment of the probability of the foundation failure.

5 Conclusions

This paper presents a 3D analysis of the bearing capacity of square and strip footings on a cohesive–frictional soil. The cohesion and friction angle were modelled using random fields based on the study of the natural soil deposit, Taranto clay. The parameters for the probability distributions and the value of the vertical SOF were adapted from existing studies. Due to the lack of information in the literature regarding the horizontal SOF, a parametric study of that value was performed for the square footing problem. For the strip footing problem, the influence of 3D modelling on the bearing capacity was analyzed and the parameter study concerned the effect of the modelled length of the foundation on the boundary value problem.

The following conclusions can be drawn from this research:

• (1) The issue of the random bearing capacity of a square or rectangular footing from the point of view of the reliability analysis constitutes a very important problem, as these kinds of foundations are very common in engineering practice. The problem requires 3D probabilistic modelling. Until recently, such a task has been beyond the computational capabilities of personal computers, and the analysis presented in this study indicates that this situation is changing.

• (2) It was shown that the modelling of a random substrate with spatially variable properties requires much more attention than in a non-random case. Apart from the obvious inability to simplify the task due to a lack of symmetry (in random modelling, the latter does not occur), the size of the domain must also be properly selected in order to avoid boundary effects for all considered realizations of random fields.

• (3) As shown for the square footing problem, the mean value for the bearing capacity increases with increasing values of horizontal SOF, but the changes are contained in a small range, and the value stabilizes at higher values of SOF.

• (4) The CoV of the square footing bearing capacity increases significantly with an increase in the horizontal SOF, but it also stabilizes for large SOF values. The results for the CoV are comparable to those obtained in a 2D analysis of the same soil (see Pieczyńska-Kozłowska et al. [30], Fig. 8).

• (5) Interesting conclusions also result from the 3D analysis of the strip foundation. The results clearly show that analyzing the strip footings in 2D overestimates the CoV of the bearing capacity, and it is probably a very conservative approach. Moreover, 3D modelling of the strip footing seems to be necessary for the precise assessment of the probability of the foundation failure.

• (6) In this analysis, almost all the features of the actual foundation were modelled. The considered footing was an embedded one, placed on a self-weighted, natural cohesive–frictional soil. The only phenomena not modelled were eccentric and horizontal loads. It seems that their introduction into the model, although associated with certain complications (i.e. modelling of additional interfaces), is not a significant problem.

The presented analysis is not without limitations. First of all, due to the relatively high computational cost, the number of simulations in each series of calculations was limited to N = 300. Although the results display some trends, this number is too small to properly estimate the probabilities of a foundation failure at their typical small values. The work associated with the modelling of the strip foundation also requires continuation (e.g. variation in the horizontal scale of fluctuation). The rapid development of computer processors as well as techniques accelerating reliability calculations provides hope that more accurate probabilistic modelling made in reasonable computation time will soon be possible.