Stochastic analyses are currently used extensively in a wide range of geotechnical applications. This is the reason of soil spatial variability caused by the geological processes that form the soil layer (e.g., Ferreira et al., 2015). The influence of soil spatial variability is examined via numerical methods based on the finite element method, finite difference method, or finite limit analysis. These methods are used to examine the impact of soil strength spatial variability on engineering structures such as foundations, slopes, or sheet pile walls (e.g., Griffiths et al., 2002; Fenton and Griffiths, 2002, 2008; Griffiths and Fenton, 2004; Srivastava et al., 2010; Kasama and Whittle, 2011; Huang et al., 2013; Simoes et al., 2014; Zhu et al, 2017; Kawa et al., 2019; Halder and Chakraborty, 2020; Pramanik et al., 2020; Żyliński et al., 2020; Kawa and Puła 2020; Viviescas et al., 2020; Pramanik et al., 2021). However, the numerical efficiency of most existing probabilistic methods is not sufficient for practical applications where three-dimensional analyses are crucial. This is true especially for three-dimensional issues for which the three-dimensional nature of soil strength spatial variability suggests that there are difficulties with its two-dimensional simplifications (e.g. plane strain conditions). This is the reason researchers investigated the possibility of providing more efficient 3-D approaches for bearing capacity analysis in the case of soil spatial variability, e.g., Chwała (2019; 2020), Li et al. (2021), Li et al., (2021). In the earlier papers by Puła and Chwała (2018) and Chwała (2019), the algorithm for using the upper bound limit theorem and Vanmarcke's spatial averaging (Vanmarcke, 1977a,b; 1983) is presented for two-dimensional and three-dimensional bearing capacity evaluation, respectively. However, the numerical efficiency achieved by Chwała (2019) can be further improved by using the so-called constant covariance matrix approach (Chwała, 2020). The constant covariance matrix approach means that the matrix is obtained using the expected values of soil strength parameters; the concept was first proposed by Puła (2004; 2007). Generally speaking, the above approaches are based on random field discretisation to a single random variable (each random variable corresponds to a specific dissipation region in a failure mechanism). These variables are correlated with the covariance matrix, and this matrix is the basis for generating the average soil strength parameters (e.g., Fenton and Griffiths, 2008; Puła and Chwała, 2015). However, the solution suggested by Chwała (2019) for three-dimensional issues is based on one iteration of the covariance matrix (for more details see the next section). The same approach is used for two-dimensional analyses by Puła and Chwała (2018). Therefore, there is still a need of verify this assumption. To study the impact of the iteration number of covariance matrix and to extend the use of the proposed method to other geotechnical applications, an iterative algorithm was developed and is described in this study. The main objective is to provide a general procedure that can be utilised for a variety of geotechnical applications, where the upper bound limit theorem (e.g., Shield and Drucker, 1953) and Vanmarcke's spatial averaging are used. Another important objective of this study is to compare the simplified approach, which is based on the constant covariance matrix with the iterative procedure that is proposed here. The iterative procedure is based on redetermining the covariance matrix in the subsequent iterations; thus, the covariance matrix corresponds to the failure geometry. As a numerical example, the proposed algorithm is used to analyse the three-dimensional bearing capacity of a shallow foundation under undrained conditions and the bearing capacity of two-layered soils. The comparison of the iterative approach with the constant matrix approach is crucial for further practical applications of the simplified algorithm for three-dimensional issues like soil sounding location optimization (Chwała, 2020; 2021). This comparison is important due to the enormous reduction in computation time. Namely, three-dimensional analyses that were performed in this study for a constant covariance matrix took about 0.5 s using a standard notebook. Therefore, there is a great potential to improve the numerical efficiency in the three-dimensional probabilistic analyses by using the constant covariance matrix approach. However, the iterative procedure is needed to examine the impact of this assumption on the resulting bearing capacity probabilistic characteristics.

The proposed procedure is devoted to random analyses of geotechnical problems in the framework of the upper bound limit theorem. Its engineering applications are directed toward foundation bearing capacity evaluation and slope stability evaluation for spatialy variable soil. The procedure is described in as general manner as possible; however, some information directed to the examples used in the study are given. Two problems are investigated employing the proposed iterative algorithm, i.e., the algorithm proposed by Chwała (2019) for 3-D bearing capacity analysis, and the algorithm proposed by Chwała and Puła (2020) and Chwała and Kawa (2021) for two-layered soil bearing capacity calculation. As mentioned in the previous section, the procedure is based on the determination of the covariance matrix, which corresponds to the actual geometry of failure. The components of the covariance matrix are established via Vanmarcke spatial averaging, and the necessary formulae have to be derived independently for each problem (those for the three-dimensional bearing capacity are given by Chwała (2019), and for two-layered soil by Chwała and Puła (2020)). As indicated earlier, the above-mentioned algorithms use one iteration of the covariance matrix, meaning that after determination of the covariance matrix (which is based on the failure geometry that was obtained for uncorrelated soil strength parameters, see details in the algorithm description below), the correlated soil strength parameters are generated. Then, they are used to find the optimal failure geometry and the corresponding final bearing capacity. However, there is an inconsistency in the above approach, which is that the optimal failure geometry differs from the failure geometry for which the covariance matrix was established. This may introduce additional uncertainty in the observed results, which can be eliminated by the iterative procedure that is proposed in this study. In the algorithm detailed below, the covariance matrix is determined as many times as the assumed number of iterations

^{2}, the type of probability density function, and the scale of fluctuation _{x}_{y}_{z}

_{2}, and the type of probability density function), an appropriate number of independent soil strength parameters _{0}^{i}_{0} are generated using a random number generator. This number is equal to the number of dissipation regions in the considered failure mechanism

_{0}_{0} (bearing capacity) is computed based on the total energy dissipation in the failure mechanism. However, it is necessary to find the minimum upper bound limit load, and the optimal failure geometry is searched using a dedicated optimisation procedure. Thus, the optimal failure geometry Δ_{0} and the corresponding upper bound limit load _{0} are determined. The optimisation procedure should be chosen individually for the considered problem, and there are no specific restrictions on this issue.

_{0}_{0}. The size of the matrix is _{1} and _{2} are the integration regions (section, surface, or volume). Based on _{0}_{0}^{C}_{1}_{0}_{0}

^{C}_{k}_{k}_{k}

_{k}_{k}. Based on _{k}_{0}^{C}_{k}_{+1}_{k}

^{C}_{k}_{+1}_{k+1} and corresponding upper bound limit load _{k}_{+1} are calculated.

_{max}_{max}

_{kmax} and _{kmax} are the optimal failure geometry and the corresponding upper bound limit load, respectively.

First, the above steps of the algorithm can be developed and adjusted in an individual pattern to be suitable for specified issues that are examined, such as slope stability or foundation bearing capacity.

Steps 2 to 10 are repeated

As discussed above, the difference between subsequent covariance matrices _{k}_{k}_{+1}_{k}_{k}_{+1} and failure geometry parameters Δ_{k,i} and Δ_{k+1,i}. The tendency can be expressed by Eq. (4).
_{max}_{1} and _{k}

As the first numerical example, the random bearing capacity of two-layered soil was analysed via the iterative algorithm described above. The optimisation procedure and formulae for the covariance matrix components were assumed according to a previous study (Chwała and Puła, 2020) and are not repeated here. Thus, the optimisation procedure that is based on the simulation annealing (Kirkpatrick et al., 1983; Kirkpatrick, 1984) was used to search for the optimal failure geometry (the geometry for which the bearing capacity reaches its minimum). The corresponding failure geometry that is adapted to the probabilistic analyses is shown in Fig. 2. Note that the failure mechanism shown in Fig. 2 is for the case of the homogenous sand layer as the top layer and spatially variable clay as a bottom layer. This scenario describes the situation of working platforms, where the top layer homogeneity results from a man-made layer. In this example, the boundary between two soil layers is assumed to be deterministic, i.e., its depth is constant in Monte Carlo analyses; however, generally, this issue is more complicated (e.g., Ghanem and Brząkała, 1996; Bagińska et al., 2020; Rainer and Szabowicz, 2020).

First, as a numerical procedure, the iterative algorithm that is denoted by ‘A’ in Fig. 1 was used. The assumed iteration number is _{k}

Eight scenarios were considered for the bearing capacity of two-layered soil. Two averaging lengths are considered (they can be interpreted as shallow foundation lengths), i.e., a _{v}_{v}_{h}_{h}_{cu} = 9 kPa and coefficients of variation _{cu} = 0.5. In this scenario, a linear trend in undrained shear strength is also considered. The included trend is 9 kPa/m starting from the bottom layer. The trend consideration is based on the recent algorithm by Chwała and Kawa (2021). For the homogenous sand layer, the friction angle

According to the algorithms presented in Fig. 1, the number of realisations _{h}_{v}

The obtained results for two-layered soil provided as bearing capacity mean values and standard deviations are shown in Fig. 4a and Fig 4b, respectively.

As shown in Fig. 4a, the bearing capacity mean values are in a narrow range, the greatest differences occur in a range for _{v}

As shown in Fig. 5, both the bearing capacity mean value and standard deviation are very close to each other for all considered values of _{h}

As the second numerical example, the random bearing capacities of square and rectangular foundations were analysed via the iterative algorithm described above. The optimisation procedure and formulae for the covariance matrix components were assumed according to a previous study (Chwała, 2019) and are not repeated here. Thus, the optimisation procedure that is based on the simulation annealing (Kirkpatrick et al., 1983; Kirkpatrick, 1984) was used to search for the optimal failure geometry (the geometry for which the bearing capacity reaches its minimum). A rough foundation base is assumed, and therefore, the corresponding failure geometry is adapted to the probabilistic analyses, and is shown in Fig. 6.

First, as a numerical procedure, similarly as in the previous example, the iterative algorithm ‘A’ (see Fig. 1) is used. The assumed iteration number is _{k}

Six scenarios were considered for the three-dimensional failure mechanism. Two shapes of the shallow foundation are assumed, i.e., square foundation (_{v}_{h}_{cu}=100 kPa and coefficients of variation _{cu}==0.5. In this scenario, a linear trend in undrained shear strength is not considered.

The number of simulations

As shown in Fig. 7, the bearing capacity mean values are in a narrow range, the greatest differences occur in a range of

As shown in Fig. 8, both the bearing capacity mean value and standard deviation are very close to each other for all considered values of

The results obtained in the two previous sections are for different applications of the random failure mechanism method (RFMM). The first concerns the bearing capacity of two-layered soil when the failure mechanism is based on two-dimensional simplification, and the second concerns the bearing capacity of the rectangular foundation in the fully three-dimensional case. In both cases, soil spatial variability is considered three-dimensional. According to the obtained results, the differences between the approach with constant covariance matrix ‘C’ (_{u}_{cu}=0.25, _{cu}=0.5 and _{cu}=1.0. A very slow increase in the differences between mean values and standard deviations of bearing capacity in both cases (_{cu} (see Fig. 9a and Fig. 9b). Despite those differances, the resulting bearing capacity coefficient of variation is almost the same for each considered _{cu} (see Fig. 9c).

All numerical results shown in this study were performed for relatively small number of Monte Carlo realisations. However, as mentioned earlier in this study the main goal is to examine the impact of using more iterations of the covariance matrix. For this particular approach the accurate determination of bearing capacity mean value and standard deviations is not necessary. Despite this, in Fig. 10 an example of such analyses is shown to demonstrate that even relatively low number of Monte Carlo realisations (

The study propose an original iterative algorithm (path ‘A’ or ‘B’ in Fig. 2) for the upper bound analyses of geotechnical problems with consideration of soil strength spatial variability. The resulting limit loads are realisations of a random variable for which probabilistic characteristics must be determined using the Monte Carlo method. The use of the iterative algorithm ensures the consistency between the failure geometry and the covariance matrix, from which the average soil strength parameters are determined - the averaging procedure is following Vanmarcke's spatial averaging (Vanmarcke, 1977a, 1977b, 1983). Moreover, the general algorithm for using a constant covariance matrix is presented in the study. Performing an analysis via the algorithm based on the constant covariance matrix (see algorithm ‘C’ in Fig. 2) provides a dramatic improvement in computational efficiency. Therefore, the possibility of using a constant covariance matrix is very promising for practical applications and three-dimensional analyses. The main objective of the study is to investigate weather these algorithms give similar results. According to the delivered examples it is shown that, the constant covariance matrix approach provides satisfactory results. However, the investigation required a comparison of the results provided by both approaches. This was done through the iterative algorithm proposed in this study. However, despite good agreement between the two considered approaches for the analysesd scenarios, i.e., bearing capacity of rectangular foundation and two-layered soil under plane strain assumption, other geotechnical applications need to be verified individually. Therefore, the conclusions are only limited to those two scenarios.

The proposed iterative algorithm was utilised in this study to evaluate the bearing capacity of shallow foundations. Two cases were considered, i.e., the bearing capacity of two-layered soil, and the three-dimensional analyses, both were performed for a variety of foundation shapes and scales of fluctuation. To summarize, the following conclusions can be drawn:

The study presents an original iterative algorithm for upper bound analyses of geotechnical problems with the inclusion of soil strength spatial variability. The proposed algorithm can be applied not only to bearing capacity evaluation, as shown in the study, but also for other applications like slope stability or retaining wall analyses. The use of the iterative procedure ensures consistency between the failure geometry and the covariance matrix; therefore, the algorithm allows the recognition and control of the impact of this issue on the resulting estimates. Basing the iterative algorithm on the upper bound theorem provides the opportunity to utilise it as a reference for other probabilistic methods.

This paper shows the possibility of applying the iterative algorithm for shallow foundation random bearing capacity evaluation. Moreover, the proposed algorithm allows the authentication of previous results (e.g., Chwała, 2019; Chwała and Puła, 2020) and indicates that for the undrained conditions, two iterations (

Scenarios analysed by the iterative algorithm were analysed again using the algorithm based on the constant covariance matrix (see the algorithm ‘C’ in Fig. 2). The obtained results are presented to allow easy comparison between both approaches (see Figs 4, 5, and Fig. 7, 8). The results provided by the algorithm, which use the constant covariance matrix concept, are very close to those obtained by the iterative procedure. This is very promising for further application of the algorithm ‘C’ in solving practical problems, e.g., optimal placement of soil soundings (Chwała, 2020; Chwała, 2021). It is also the most important result of the study that opens up the possibility of using the constant covariance matrix. Thus, a dramatic improvement in the computational efficiency is possible (analysing one realisation of a three-dimensional bearing capacity problem including soil strength spatial variability, takes about 0.5 s for a standard notebook). To summarize, the algorithm based on the constant covariance matrix can be used for three-dimensional random bearing capacity problems and other issues that require high computational efficiency.

^{th} ICSMGE; Paris 2013;731–734.

Effect of Specimens’ Height to Diameter Ratio on Unconfined Compressive Strength of Cohesive Soil Effect of fibre content on the geotechnical properties of peat The new railway hybrid bridge in Dąbrowa Górnicza: innovative concept using new design method and results of load tests The Temperature Field Effect on Dynamic Stability Response of Three-layered Annular Plates for Different Ratios of Imperfection Safety of Steel Arch Support Operation During Rock Bursts Under Explosive Atmosphere Conditions