Evaluation of shallow foundation bearing capacity in the case of a two-layered soil and spatial variability in soil strength parameters

In this study, a probabilistic approach for evaluating the bearing capacity of surface footings is discussed. The evaluation is based on a kinematic approach. The considered substrate consists of two different layers of soil: a top layer formed of medium or dense sand followed by a layer of soft clay. The sand layer is assumed to be homogenous, whereas the undrained shear strength of the soft clay layer is assumed to be spatially variable, described by a lognormal random field. The random field is discretized according to Vanmarcke’s spatial averaging along dissipation regions in the considered failure mechanism. The mechanism utilizes plane strain conditions; however, due to consideration of the soil spatial variability in three dimensions, the impact of the length of the foundation on the random bearing capacity evaluation is considered in this study. As a result of the discretization procedure, a set of correlated random variables is obtained (each associated with an individual dissipation region in the failure mechanism). A series of numerical analyses are performed for two thicknesses of the first layer and a set of anisotropic correlation structures for the spatial variability of the undrained shear strength. The proposed method is computationally efficient and allows consideration of three-dimensional spatial variability in soil strength properties. The results are discussed and compared with those obtained by other methods.


Introduction
Bearing capacity evaluations of shallow foundations in the case of spatially variable soil are currently used primarily for two-dimensional analysis. One of the reasons for this situation is numerical efficiency, which limits the use of three-dimensional finite element formulations for random analyses. However, recent attempts have been made to solve such issues, e.g., [1]. Moreover, most existing probabilistic methods have been applied to single-layered soil. However, from an engineering perspective, two-layered soil is also a highly probable situation. An especially important situation is when the bottom layer is a soft soil and the top layer is a strong soil. This scenario has been extensively examined in the case of deterministic analyses (e.g., [2][3][4][5]); however, there have been very few studies on random analyses of two-layered soil (e.g., [6]). It is worth mentioning that it is possible to perform these types of two-dimensional analyses by using OptumCE software. a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 However, the lack of random approaches that are efficient and capable of introducing three-dimensional analyses was the motivation to propose an original approach that utilizes a kinematic approach in conjunction with Vanmarcke spatial averaging [7][8][9]. An analogous algorithm for plane strain conditions and a multiblock failure mechanism for a one-layered soil was proposed by Puła and Chwała [10]. The present study is also motivated by the engineering application to the problem of the bearing capacity of working platforms, which are usually constructed over soft subsoil to ensure safe traffic of heavy trucks. A working platform usually consists of compacted cohesionless material. Due to the controlled compaction procedure, the relatively homogenous soil layer is obtained. However, the spatial variability of the natural soil layer situated below the man-made layer can play a significant role in bearing capacity estimation. This situation is reflected in the numerical analyses. The importance of the problem can also be illustrated by investigating failures of working platforms, e.g., Channel Tunnel Rail Link Contract 310 [11] or Lublin (eastern Poland) road bypass (e.g., [12]).
In this study, as a deterministic background, a failure mechanism proposed by Michałowski and Shi [3] is assumed. The selected failure mechanism is dedicated for two-layered soil, i.e., the top layer is medium or dense sand, and the bottom layer is soft clay. According to the above described assumption, the spatial variability in soil strength is considered only in the bottom layer; thus, this approach characterizes the spatial variability in the undrained shear strength. The problem of designing working platforms is coded in BR 470 [13]. This study provides an efficient method for analysing such engineering problems and demonstrates that the random bearing capacity of two-layered soil in the case of spatial variability in the bottom layer can be efficiently evaluated by the developed kinematic approach. Moreover, despite the use of a two-dimensional mechanism, the spatial variability of the undrained shear strength in the soft clay layer is considered in three dimensions. As a result, an assessment of the impact of foundation length on random bearing capacity can be estimated. The importance of performing three-dimensional analyses is noted in this study. This is especially true if soil spatial variability is considered; thus, the spatial variability has a three-dimensional structure, and using two-dimensional simplifications can significantly influence the final results. A numerical algorithm is proposed in this study. According to this numerical procedure, numerous two-layered soil scenarios are analysed. Moreover, a simple approach to separate the impact of the thickness of the homogenous layer on random bearing capacity estimates from the impact of the spatial variability of the bottom layer and foundation length is provided in the study. The obtained results are shown and discussed in section 5. In summary, the objective of this study was to propose an efficient approach for threedimensional bearing capacity analysis of two-layered soil. The considered situation reflects the working platform scenario in which the top layer is a man-made layer of compacted cohesionless soil (in this study assumed to be homogenous) and the bottom layer is a natural, spatially variable, cohesive soil.

Failure mechanism description for random analyses
The deterministic failure mechanism proposed by Michałowski and Shi [3] is used in this study. The geometry of the failure mechanism was modified to be suitable for random analyses. As a result, a new formula for bearing capacity is introduced. A detailed description of the random failure mechanism is provided below. The first assumption concerns the homogeneity of the top soil layer (note that in this study, a sand layer is considered to be the top layer). In the case of a working platform, the upper layer is practically a man-made material. Therefore, the spatial variability of this layer is neglected, and only a deterministic case is provided. As a result, the spatial variability in soil properties is considered for the bottom layer (purely cohesive soil). Moreover, for the bearing capacity estimations performed in this study, the bottom layer is assumed to be weightless. In the case of undrained conditions, the weight of the soil can be neglected. Moreover, as shown in Chwała [14] and also in the case of random analyses, the consideration of soil weight has a very limited impact on the final bearing capacity estimations. As in the deterministic version of the failure mechanism [3], the random failure mechanism proposed in this study consists of rigid blocks of soil that move with respect to each other. This failure mechanism is constructed in a particular manner so that there is no velocity discontinuity along the line that separates the top soil layer from the bottom layer. For the numerical analyses performed in this study, the failure mechanism of the 11th block is assumed. An example geometry and the rigid block numbering convention are shown in Fig 1. Note that the block shown in Fig 1 is a single rigid block, i.e., the block denoted by O 1 P 1 A 1 A 2 P 2 O 1 defines a rigid block. Due to assumption of the homogeneity of the sand layer (top layer), all velocity jump vectors within this layer are inclined to the slip lines at the same angle that is equal to the assumed internal friction angle of sand φ (this angle is shown in Fig  2). Note that the velocities shown in The first noticeable difference when comparing the failure mechanism shown in Fig 2 with the deterministic failure mechanism is that in the random case, the failure geometry is no longer symmetric. Due to the different soil strength properties for each velocity discontinuity, the failure geometry is asymmetrical with respect to the vertical line passing through the foundation centre. An example random failure geometry is shown in detail in Fig 2. Based on Fig 2  and the upper bound theorem, the formula for bearing capacity can be derived. As an illustrative example, some parts of the bearing capacity formula are derived below. Based on the upper bound theorem, the upper limit of the bearing capacity can be found by comparison with the dissipation energy along velocity discontinuities and gravitational forces. Note that no overburden pressure is considered (a foundation is placed on the soil surface, as shown in Figs 1 and 2). According to results presented by Chwała [14], the self-weight of soil in the case of purely cohesive soil can be neglected when the spatial variability of soil strength is considered. Rigid block numbering convention for the random version of the failure mechanism for two-layered soil. Note that t is the thickness of the top layer (homogenous sand) and that the bottom layer is assumed to be spatially variable clay. https://doi.org/10.1371/journal.pone.0231992.g001

PLOS ONE
As an example, the energy dissipated along the velocity discontinuity line P 1 A 1 can be calculated as shown in Eq (1).
where c 1 and l 1 are the cohesion and dissipation line length, respectively, and v 1 is the velocity jump vector. All three values are shown in Fig 2. By using an analogous method, the remaining formulas for energy dissipation within the bottom layer can be derived. Moreover, because of the cohesionless character of the top layer, there is no energy dissipation within it (c = 0). Therefore, in the considered failure mechanism, the energy dissipates along sixteen lines within the bottom layer (see Fig 2). Another factor that influences the final bearing capacity estimation is the power of the gravitational forces on the top layer. The formula for region O 1 P 5 C 1 is demonstrated in Eq (2).
where F 6 is the area of the rigid block indexed by 6 (see Fig 1). Note that in the case of block 5, the corresponding area F 5 is the area of block 5 within the top layer. The velocity vector in brackets [v 10 ] denotes the vertical component of the velocity jump vector v 10 . The unit weight of the sand layer is denoted by γ. By using an analogous method, the remaining formulas for the power of gravitational forces within the top layer can be derived. The power of the limit load that the foundation can be subjected to can be calculated by multiplying the bearing capacity p by the velocity v 0 of the rigid block beneath the foundation (the foundation moves with the same velocity that is vertically oriented). Note that the formulas provided above assume by default that the calculations are performed for a = 1.0 m (see Fig 2); however, if other foundation lengths are considered, the corresponding value of a must be included. By summing all terms that relate to energy dissipation and gravitational forces and by transforming the resulting equation due to p, the following formula can be obtained: where, as indicated earlier, v i is the magnitude of the velocity jump vector along the slip line i, [v i ] denotes the vertical component of velocity v i (negative if upward), and l i is the length of the slip surface i. In Eq (3), the undrained shear strengths are denoted simply by c i , where the index i is the same as that for the slip line lengths l i (see Fig 2). Design charts for calculating bearing capacity were provided in a paper by Michałowski and Shi [3]. As shown in Fig 2, when increasing the depth of the bottom layer (t), the critical value of t can be found. This critical value means that the failure mechanism is no longer deep (spread out in both soil layers) but becomes much shallower (spread out only in sand). This issue was described in detail by Michałowski and Shi [3], who showed that the critical depth of the weak layer depends on the angle of internal friction φ of the sand and the ratio c u /γb. Nevertheless, this issue is not considered in this study; thus, the examined scenarios relate to the case in which a deep failure mechanism occurs (as shown in Figs 1 or 2). The investigated situation is adequate for working platforms in which the top sand layer is relatively thin.
The formula shown in Eq (3) is crucial for determining the bearing capacity in the case of spatial variability in the strength properties of the bottom layer by applying different undrained shear strengths on each dissipation region. However, it is essential to properly calculate these undrained shear strengths. The method used in this study for that purpose is detailed in the next section. Moreover, Eq (3) provides the bearing capacity if the failure geometry is known (see Fig 2). Therefore, an optimization procedure is necessary to find the minimum value of the limit load because the method is based on the upper bound approach.

Spatial averaging
As indicated in the previous section, the method for determining average undrained shear strengths is a crucial element of the proposed approach. Because the sand layer is assumed to be non-random, the spatial averaging is subject only to the random field X that describes the spatial variability of the bottom layer. For this purpose, the method proposed by Puła and Chwała [15] is used here. Note that the approach described in the mentioned paper concerns a one-layered soil and Prandtl failure mechanism. However, with some modifications and adaptations, this approach can be utilized for the scenario considered in this study. The modified version of the spatial averaging procedure is described below. Additionally, to maintain the clarity of this paper, some general information is also provided below.
The approach is based on Vanmarcke's spatial averaging concept [7][8][9]. Generally, the procedure replaces the need to generate a random field X by generating a set of correlated random variables (random vector X V ). Those random variables (X V,1 , X V,2 , . . ., X V,n ) or random vector components correspond to the dissipation regions resulting from the considered failure mechanism. Due to the assumed stationarity of the random field, after the averaging procedure, the mean values of the random vector components are preserved; however, their variances are subject to reduction. The average is calculated as shown in Eq (4): where V i is the averaging domain that corresponds to the dissipation region and X is the initial random field. Because spatial averaging is performed in three dimensions, the dissipation regions in this study are rectangles (see Fig 2). The variance of the new random variable X Vi is given by: where s 2 X is the variance of the random field X. According to earlier experience (e.g., [16,17]), in this study, the field X is assumed to be a lognormal random field. However, both isotropic and anisotropic correlation structures are examined. As a correlation function within the random field, the Gaussian covariance function was selected and used in the numerical analyses, as shown in Eq (6).
where θ x , θ y and θ z are fluctuation scales [7,17]. To obtain the covariance function from Eq (5), the right-hand side of the equation must be multiplied by s 2 X , i.e., RðDx; Dy; DzÞ ¼ s 2 X rðDx; Dy; DzÞ. Note that in this study, (5) is applied to the three-dimensional scenario. Despite the assumed plane strain conditions, spatial averaging is performed in three dimensions. The coordinate system is shown in Fig 2. However, in further analyses, both horizontal fluctuation scales are assumed to be equal. Moreover, they are not distinguished and are denoted as θ h (θ h = θ x = θ y ). The fluctuation scale along the z direction is called the vertical fluctuation scale θ v (θ v = θ z ). According to the failure mechanism shown in Fig 2, there are 16 regions in which spatial averaging is necessary. To determine these regions, a covariance matrix that describes the mutual correlations between each pair of regions is needed. To calculate the components of the covariance matrix, dedicated formulas must be derived. Derivation is possible due to the equations shown in Eq (7) (see [15,18]).
The above equation determines the covariance between two random variables X Vi and X Vj resulting from the averaging (Eq (4)) corresponding to the dissipation regions V i and V j (here, rectangles are shown in Fig 2; e.g., P 1 A 1 P' 1 A' 1 ). Note that i, j = 1, . . ., 16. As a result, the final size of the covariance matrix is 16×16. The covariance matrix is symmetric and positively definite. Below, some formulas for the covariance matrix components are derived. First, consider the dissipation region P 1 A 1 P' 1 A' 1 (points P' 1 and A' 1 are not shown in Fig 2, but their coordinates are the same as those of P 1 and A 1 , respectively, except that the y coordinates are shifted by α). Fig 2 shows the angle α (coloured red), which is used for parametric representation of the region P 1 A 1 P' 1 A' 1 . Indeed, the considered region can be expressed by using the following parametrization (Eq (8)): Eq.d, by using the following parametrization the considered region can be expresses to be derived. is needed.
By substituting Eqs (8) and (6) into Eq (7), the following formula for the new variance of region P 1 A 1 P' 1 A' 1 is obtained: Note that in Eq (9), the region P 1 A 1 P' 1 A' 1 is denoted more simply as P 1 A 1 ; this convention is also used later. Other variances of regions described by points P i and A i can be derived analogously.
In the case of region A 1 A 2 , the transformation shown in Eq (10) can be utilized as follows: where the angle δ is shown in Fig 2 (coloured red).
Note that in Eqs (9) and (11), the coordinates along each axis are subtracted from each other. As a result, the coordinates of points P 1 and A 1 are reduced. However, the situation is different when the covariance between two different regions is calculated. To illustrate this situation, the formula for the covariance between regions P 1 A 1 and A 1 A 2 is given in Eq (12).
Based on an analogous procedure, the other covariance matrix elements were derived. To enable presentation of the covariance matrix, a simplification in the naming convention is introduced. First, start from the components considered above, i.e., w 1 ¼ Namely, the variances are denoted by w i , where i is the index of the considered region (see Fig 2); the covariances are denoted by k ij , where i 6 ¼ j. As a result, the covariance matrix is given in Eq (13).
In the presented method, the covariance matrix size depends on the number of dissipation regions resulting from the failure mechanism. The size given in Eq (13) is valid for the failure mechanism of the 11th block shown in Fig 2. However, if a more rigid block is considered, the size of the covariance matrix is correspondingly larger.
The method for determining the covariance matrix for the considered failure mechanism is provided above. Next, based on C X , the average set of undrained shear strengths is computed by the method described in Appendix A based on the approach proposed by Puła and Chwała [15]. Through the algorithm, the sixteen independent values of undrained shear strengths (c u,1 , . . ., C u , 16 ), each dedicated to a specific dissipation region, are transformed to the correlated values (c u;1 ; . . . ; c u; 16 ). The transformation is based on the Cholesky decomposition of the covariance matrix (Eq (13)). More details on this issue are also provided in the numerical algorithm section.

Simulated annealing optimization scheme
As mentioned earlier, the procedure for calculating bearing capacity described in section 2.1 is applicable for a known failure geometry. However, before the covariance matrix is determined, the optimal failure geometry must be found. Here, optimal means the failure geometry that provides the lowest possible bearing capacity for the specified soil strength parameters. Therefore, the optimization procedure is crucial to use the multiblock failure mechanism shown in Fig 2. As a basic optimization procedure, a simulated annealing scheme is utilized in this study [19,20]. Moreover, earlier authors' experiences with the application of simulated annealing to optimize the geometry of kinematic failure mechanisms noted that this approach could be used here [10,21]. As a result, the approach described below is created. The flow chart of the procedure is shown in Fig 3; to facilitate navigation, the step numbers have been added in the below description. Generally, the optimization method is based on slight changes in the failure geometry in subsequent simulations. The method starts with the initial geometry parameters (step 1) from which the first bearing capacity is calculated (step 2). After that, the control parameter values are set (step 3). For each simulation, the corresponding bearing capacity is calculated using Eq (3). One bearing capacity is called the current value p c , and this value is compared with those from the subsequent simulations p n . This method allows for accepting a worse solution than the current one (p c < p n ); however, the probability of such a scenario P a decreases during the simulation process, and at the end of the simulation, this probability is nearly zero (only better solutions are accepted). Note that here, the better solution means the failure geometry that provides lower bearing capacity estimation. The concept of accepting worse solutions is due to overcoming local extrema; therefore, the final result is an approximation of the global extremum. In the simulated annealing scheme, there are some controlling parameters that control the simulation process (step 3). The formula used to calculate the acceptance probability is shown in Eq (14).
The parameter T is responsible for decreasing the acceptance probability P a during the simulation process. Therefore, T also decreases during simulation. Note that Eq (14) is used if p c < p n ; however, if p n < p c , the formula given in Eq (14) reaches a value greater than one, and as a result, a better solution is always accepted (step 8). According to the literature (e.g., [22]), to effectively set the optimization method, the average initial value of P a should be approximately 0.5. The crucial element of the described optimization procedure is a procedure for generating a slightly different failure geometry based on the current one. This is a so-called neighbouring set of geometric parameters (step 6). This procedure is used to generate each failure geometry considered by the simulation process (omitting the first geometry that is determined arbitrarily). The geometry of the failure mechanism shown in Fig 2 has 18 degrees of freedom. Therefore, those 18 parameters are the arguments of the objective function that is intended to minimize Eq (3). As a result, the arguments of the procedure for determining a neighbouring failure geometry are those 18 parameters. The following parameters are selected for this purpose: x C 1 and x C 2 coordinates; x P i coordinates, i = 1, . . ., 10; and lengths l 2 , l 3 , l 4 , l 10 , l 11 and l 12 . Each geometric parameter is changed in a random order by increasing its initial value by an amount generated from a uniform distribution U[−0.01 m, 0.01 m]. Moreover, the procedure takes into account additional limitations on the values of the changed parameters, for example, to prevent overlapping of adjacent blocks, namely, to ensure that, e.g., x P 7 > x P 6 (see Fig 2). As mentioned earlier, for a new set of geometric parameters, the bearing capacity p n is calculated using Eq (3) (step 7) and compared with the current value p c (step 8). The acceptance probability of a worse solution is modified by the parameter T; however, some simulations are performed under constant P a . In this study, for each P a level, z = 2000 simulations are carried out (step 5). Then, P a is decreased by modifying the value of T, namely, by multiplying T by an additional parameter: αT (step 9). The parameter α is kept constant during the simulation process and based on the literature and earlier experiences is assumed to be equal to 0.9. The simulation process continues as long as T > T min (step 4). In this study, T min = 10 −8 . When T reaches T min , the algorithm ends, and the current bearing capacity p c is treated as the optimal bearing capacity (step 10). Note that the values of the parameters that control the simulation process have to be selected individually for the considered issue.
According to the performed tests, the parameters chosen for the purpose of this study are definitely appropriate in terms of the achieved accuracy. The geometry of the failure mechanism obtained by using the optimization procedure described above for a non-random case was identical to the geometry obtained by Michałowski and Shi [3]. One example of this comparison is shown in Fig 4. For the initial geometry, the bearing capacity is approximately 210 kN/m. However, at the beginning stage of the simulation process, the value reaches over 280 kN/m. Despite this, during a simulation process with decreasing P a , the current bearing capacity p c moves towards the optimal value.
The results obtained for the optimization procedure were also compared with the results obtained using finite element limit analysis implemented in OptumG2 [23] software for different numbers of finite elements and calculation methods. The obtained results are juxtaposed in Table 1. The results were obtained under plane strain conditions. To enable their comparison to the three-dimensional case, additional deterministic analyses were carried out by using ZSoil software [24]. In the case of t = 0.

Algorithm summary
In sections 2.1, 2.2 and 2.3, all the necessary components to build the numerical procedure are described. The connection of those elements is provided in the next section to find the probabilistic characteristic of the random bearing capacity. As a governing framework of the numerical procedure, a Monte Carlo method is used. Moreover, the flow chart of the algorithm layout is provided in the next section.

Numerical procedure
The numerical procedure for a one Monte Carlo simulation is provided below (see also Fig 5): Step I introduces the values of the initial parameters. Set i = 1.
Step II starts Monte Carlo loop. If i � N follow the steps below. If i > N go to step VII.
Step III generates 16 independent values of undrained shear strengths via a pseudo-random number generator Y c u ¼ ðc Y;u;1 ; . . . ; c Y;u;16 ). All values are generated from the underlying normal distribution (m Y;c u and s Y;c u ) for the initial stationary lognormal random field (m c u and s c u ).
Step IV finds the optimal bearing capacity based on the soil parameters from step III transformed to the lognormal distribution X c u ¼ ðc X;u;1 ; . . . ; c X;u;16 Þ by using the optimization procedure described in section 2.3.
Step V determines the correlated undrained shear strengths by transforming the independent values of undrained shear strength obtained in step III to correlated values using the covariance matrix (Eq (13)) determined for the optimal failure geometry from step IV. For this purpose, the algorithm shown in Appendix A is used. This algorithm is based on the Cholesky decomposition (e.g., [25]) of the covariance matrix C X . Generally, the correlated parameters are obtained by calculating the product of a vector of uncorrelated parameters obtained in step III that was standardized and a triangular matrix that is the result of the Cholesky decomposition of C X . Finally, a set of correlated values for the undrained shear strength is obtained (c u;1 ; . . . ; c u; 16 ). These values are correlated by the covariance matrix C X .
Step VI calculates the final bearing capacity by using the correlated undrained shear strengths (c u;1 ; . . . ; c u;16 ) obtained in Step V. For this purpose, the optimization procedure from section 2.3 is used again. The resulting failure geometry and the corresponding bearing capacity are returned as the final results of the algorithm. Set i = i + 1, and go to step II.
Step VII provides the set of N bearing capacity estimates.  estimations for b = 1.0 m, γ 18 kN/m 3 , φ = 35˚and two thicknesses of the sand layer: t = 0.5 m and t = 1.0 m. In this study, for all performed analyses, the number of Monte Carlo simulations N is equal to 1000. The proposed approach is computationally efficient; namely, one three-dimensional simulation takes approximately 1 minute for one core of a standard notebook.

PLOS ONE
The numerical procedure detailed above allows consideration of three-dimensional soil spatial variability in the bearing capacity estimations for two-layered soil. Unfortunately, there are no published results that can be compared with the approach proposed in this study. However, this is possible in the case of plane strain conditions. In this particular case, the random finite limit analysis can be used for a comparison by using OptumG2 software. The comparison is made for m c u ¼ 20 kPa and v c u ¼ 0:25 for two thicknesses of the top layer: t = 0.5 m and t = 1.0 m. The results for θ v = 1.0 m plotted as a function of θ h /θ v are shown in Fig 6. The observed differences between both approaches are acceptable. Moreover, to some extent, the differences can be explained by the different correlation function assumed in the analyses performed by OptumG2 software (the Markovian one).

Considered scenarios
The main objective of this study is to assess the impact of a homogenous sand layer (top layer) on the overall random bearing capacity estimation in the case of a two-layered soil. A spatially variable, purely cohesive soil was assumed as the bottom layer. Moreover, the considered scenarios were selected to ensure creation of the failure mechanism shown in Fig 2. Therefore, relatively small thicknesses of the sand layer were examined. These thicknesses are suitable for working platform problems. The impact of vertical and horizontal fluctuation scales on random bearing capacity was investigated. Three vertical fluctuation scales were selected for numerical analyses, namely, θ v = 0.25 m, θ v = 0.50 m and θ v = 1.00 m. These values were recently indicated by many authors, e.g., [26]. In the case of horizontal fluctuation scales, the values of θ h = θ v , θ h = 5θ v , θ h = 10θ v and θ h = 20θ v were examined. The two-dimensional failure mechanism was extended to three dimensions (Fig 2) to enable consideration of three- Numerical analyses of all of the abovementioned parameter combinations were performed for two thicknesses of the top layer, i.e., t = 0.5 m and t = 1.0 m. The information for all scenarios that considered the impact of fluctuation scales is juxtaposed in Table 2.
Moreover, for comparison purposes, scenarios with assumed infinite horizontal and vertical fluctuation scales were examined. By assuming θ v = θ h = 1, a simple situation where the undrained shear strength is described by a single random variable is obtained. In such a case, there is no spatial averaging; therefore, the obtained standard deviations or coefficients of variation of bearing capacity are conservative estimates of the true values. These scenarios were intended to be analysed to separate the influence of dumping the spatial variability impact on the final bearing capacity coefficient of variation for the deterministic top layer of sand.

Results
Based on the Monte Carlo analyses, the mean values of the resulting bearing capacity μ p and corresponding bearing capacity coefficients of variation were obtained. In Fig 7,  This effect can be a result of using a lognormal probability distribution of c u ; a similar effect was described in [27]. In Fig 7a, 7b and  In each case, a longer averaging domain provides a higher μ p . This can also be explained by the abovementioned effect of obtaining a lower μ p for a greater coefficient of variation c u . Thus, a longer averaging domain a provides greater averaging; therefore, the overall coefficient of variation c u decreases as the variances defined in Eq (5) are more significantly reduced. It is worth mentioning that all obtained bearing capacity mean values are below the values of the deterministic scenario that is denoted by the horizontal dashed line. Note that the relative differences in the obtained μ p values are not very significant (see values on the vertical axis). Analogous observations can be made for t = 1.0 m (Fig 8). However, the relative differences in the obtained μ p values are smaller than those for t = 0.5 m. This demonstrates the effect of damping the impact of soil spatial variability on bearing capacity with increasing thickness of the homogenous layer t. This effect is described later in more detail. Very significant differences among the considered scenarios are observed in the coefficient of variation of bearing capacity v p . The results for v c u ¼ 0:25 are shown in Fig 9. The horizontal dashed and dotted lines shown in Fig 6 indicate v p for special scenarios. The dashed line indicates the coefficient of variation of bearing capacity for t = 0 and infinite values of the vertical and horizontal fluctuation scales. In this case, the considered issue is degenerated into a one-layered purely cohesive soil with no spatial variability and v c u ¼ 0:25. Therefore, the resulting coefficient of variation of bearing capacity v p ¼ v c u ¼ 0:25 (as shown in Fig 9). However, the dotted lines in Fig 9a and 9b

PLOS ONE
further decrease is observed for t = 1.0 m, for which v p,t 0.157. Thus, it is easy to understand why the v p values observed for scenarios with finite fluctuation scales (see Table 1) are lower than v p,t . As shown in Fig 9, a shorter vertical fluctuation scale provides a smaller value of v p . The dotted lines represent the limit values of v p when the fluctuation scales reach infinity. The information is intentionally presented in the form shown in Fig 9 to separate the impact of the homogenous sand layer from the impact of fluctuation scales on the resulting bearing capacity coefficients of variation. As expected, for a longer averaging domain a, a lower v p is obtained. The observed differences between the extreme values, i.e., a = 1 m and a = 10 m, are significant, and they indicate the necessity of considering the real dimension of the structure if soil spatial variability is taken into account. In such situations, the two-dimensional simplification of three-dimensional issues provides more conservative bearing capacity estimates. This conclusion has recently been noted [1,21,29].
Results analogous to those shown in Fig 9 for v c u ¼ 0:5 are shown in Fig 10. The horizontal dashed and dotted lines have the same meaning as those in Fig 9. Both Figs 9 and 10 indicate that calculating v p,t can be a very easy method to conservatively estimate the random bearing capacity in the case of a lack of information about fluctuation scales. The reduction in the bearing capacity coefficient of variation is a result of considering only the homogenous soil layer. Therefore, it can easily be estimated. To illustrate that possibility, Figs 11 and 12 show the results obtained by analysing scenarios with θ v = θ h = 1. First, in Fig 11, the obtained bearing capacity mean values are plotted as a function of the homogenous sand layer thickness t.  can be used to estimate the reduction factor in the bearing capacity coefficient of variation that results from the thickness of the homogenous sand layer.

Summary and conclusions
An original adaptation of the two-dimensional failure mechanism for two-layered soil proposed by Michałowski and Shi [3] to random analyses is shown in the present study. The adaptation was necessary for numerical analyses that are based on the upper bound approach and are dedicated to the scenario of a spatially variable bottom layer. A method of Vanmarcke's spatial averaging in conjunction with a kinematic failure mechanism [15] was used in this study. The use of this method requires its adaptation to the considered scenario. As a result, an original approach for random analysis of the bearing capacity of two-layered soil and threedimensional spatial averaging was proposed in this study. The created algorithm is capable of considering all important parameters, such as geometrical parameters, thickness of the top layer, its unit density and fluctuation scales, that describe the variability of undrained shear strength along the vertical and horizontal directions. The proposed approach is very efficient (one three-dimensional simulation takes approximately 1 minute for one core of a standard notebook). The method is dedicated to a situation in which the top layer is homogenous sand (generally cohesionless soil) and the bottom layer is purely cohesive soil that is spatially variable. This arrangement of the soil layers has practical meaning because it reflects the conditions in which working platforms operate, where the top layer is man-made, and the assumption that it is homogenous is justified. The proposed approach is in accordance with the upper bound theorem, and thus, the obtained bearing capacity must be treated as the upper limit of the true limit load. The proposed approach is limited to bearing capacity estimations and is not dedicated to serviceability limit states, e.g., Johari et al. [30]. In the performed analyses, seismic loadings are not examined, and thus, the obtained results are reliable for a static scenario (this can influence the obtained results, e.g., [31]). However, this type of analyses will be explored in future research.
Based on the results presented in section 5 and the description of the proposed algorithm given in section 2 and section 3, the following conclusions can be drawn: 1. The algorithm proposed in this study demonstrates that the method of joining Vanmarcke's spatial averaging with kinematic failure mechanisms can deliver an efficient approach for analysing two-layered soil problems if the bottom layer is spatially variable. Further efficiency improvement is possible if a constant covariance matrix is applied [32]. 2. The obtained results indicate that the observed reduction in the bearing capacity coefficient of variation (see Figs 6 and 7) is caused mostly by three factors: thickness of the homogenous sand layer, length of the averaging domain a, and fluctuation scales that characterize the spatial variability of the undrained shear strength within the bottom layer. The impact of surcharge load is not considered in the performed analyses. The soft layer is assumed to be weightless (see arguments given in section 2.1); however, in future studies, it would not be very difficult to extend the analyses by soil weight consideration.
3. Three-dimensional issues are of special importance when soil spatial variability is taken into account in the numerical analyses. Two-dimensional simplifications provide a more conservative approach that can be significant in some scenarios, as shown in Figs 6 and 7.
4. As shown in the results section, the impact of dumping the feature of the thickness of the homogenous sand layer can be separated from the other impacts detailed in the second conclusion. An example of this separation is shown in Fig 9, where the reduction in the bearing capacity coefficient of variation is caused only by the homogenous sand layer. As expected, the reduction is greater for a greater sand thickness. This simple approach can be utilized in practical applications for creating graphs that allow nearly instant determination of the reduction in the bearing capacity coefficient of variation.
The obtained components of T c u ¼ ðc u;1 ; . . . ; c u;16 Þ are used in step IV in the algorithm given in section 3.