Random failure mechanism method for assessment of working platform bearing capacity with a linear trend in undrained shear strength

Marcin Chwa?a, Marek Kawa

Department of Civil Engineering, Wroclaw University of Science and Technology, Wybrzez˙e Wyspian′skiego 27, Wroc?aw, 50-370, Poland

Keywords:Two-layered soil Random bearing capacity Fluctuation scale Kinematical approach Upper bound Random finite difference method (RFDM)Working platform

ABSTRACT A bearing capacity evaluation for the surface strip foundation on a working platform modelled on a twolayered substrate is considered in the study. The upper layer is assumed as man-made and wellcontrolled and thus non-variable. The lower layer modelling natural cohesive soil is subjected to spatial variability of undrained shear strength.The random failure mechanism method(RFMM)is used to evaluate the bearing capacity. This approach employs a kinematic assessment of the critical load and incorporates the averaging of three-dimensional(3D)random field along dissipation surfaces that result from the failure mechanism geometry. A novel version of the approach considering an additional linear trend of undrained shear strength in the spatially variable layer is proposed. The high efficiency of the RFMM algorithm is preserved. The influences of foundation length, trend slope in the spatially variable layer, fluctuation scales, and thickness of the homogenous sand layer on the resulting bearing capacity evaluations are analysed. Moreover, for selected cases, verification of the RFMM based assessment obtained using random finite difference method (RFDM) based on 3D analysis is provided. Two types of analyses are performed using RFDM based on associated and non-associated flow rules. For associated flow rule which corresponds to RFMM, the RFMM is conservative and efficient and thus it seems preferable. However, if RFDM employs non-associated flow rule (much lower dilation angle for sand layer),the efficient RFMM is no longer conservative. For this situation, a combined approach that improves the efficiency of the numerical method is suggested.

1. Introduction

Probabilistic analyses are currently being extensively investigated fields of geotechnical engineering. This is mainly associated with the need for consideration of the risk corresponding to soil spatial variability.This need is directly indicated in the standard ISO 2394:2015(2015).The currently popular approach is to model this variability by identifying an autocorrelation function and then to use this function in modelling the soil properties by random fields(Vanmarcke, 2010). With such an approach, the probability of failure of a geotechnical structure can be assessed within a Monte Carlo framework, by generating an appropriate number of realisations of a random field; for each of these realisations solving the problem and then analysing the distributions of the results obtained. The random finite element method (RFEM) proposed by Fenton and Griffiths (Griffiths and Fenton, 2001; Fenton and Griffiths, 2003) is a particularly often employed representation of this approach.Although some other approaches,which incorporate response surface method or recently machine learning, are also used (e.g. Wang et al., 2020), the RFEM is still a very popular and universal technique.It has been applied to the probabilistic analysis of various geotechnical problems such as bearing capacity and settlement of shallow foundations(e.g.Fenton and Griffiths,2003;Pieczy′nska-Koz?owska et al., 2015; Johari and Sabzi, 2017), slope stability (e.g. Griffiths and Fenton, 2004; Jha and Ching, 2013a),retaining wall deflection(e.g.Luo et al.,2018;Kawa et al.,2021),as well as more specific problems like stochastic analysis of internal forces in piled-raft foundations (e.g. Johari and Talebi, 2021) or stability of tunnelling perpendicularly beneath an existing tunnel(e.g. Chen et al., 2019). Although the RFEM approach has many advantages (among others, it can be used to calculate the probability of failure of any numerically modelled geotechnical structure),its main disadvantage is the relatively long computation time resulting directly from the computation time spent on individual realisation.This is particularly important from the point of view of failure probability estimation - small probabilities of failure that are demanded by the geotechnical standards require a large number of realisations(determined e.g.using the Chebyshev inequality)to be precisely assessed. Thus, the accurate determination of the probability of foundation failure at an acceptable level requires a large amount of computation time,which makes the RFEM method difficult to apply in engineering practice.

In the case of bearing capacity evaluation for strip foundation,the plane-strain state and stationary random fields are usually assumed(Fenton and Griffiths,2003;Pieczy′nska-Koz?owska et al.,2015).This assumption was also used in the earlier works of the authors (Pu?a and Chwa?a, 2015, 2018; Chwa?a, 2019). The random failure mechanism method(RFMM)developed in these works,which is based on the limit state upper bound theorem,allowed to significantly reduce the computation time (compared to RFEM) necessary to assess the probability of foundation failure.

It should be noted that while the assumption of plane-strain condition is correct for deterministic analysis, in the case of employing random fields, it is a large simplification. This is due to soil spatial variability being three-dimensional (3D) phenomena.The change of properties in the direction along the foundation,which can significantly influence the results,cannot be reproduced using a two-dimensional(2D)random field representation.For this reason (as shown e.g. by Kawa and Pu?a, 2020), in the case of spatially variable soil,3D analyses are essential for strip foundation bearing capacity evaluation. On the other hand, 3D analyses are particularly time-consuming for a finite element or finite difference based approach; therefore, there are still only a few studies that apply 3D analyses for considering soil spatial variability. Such analysis was already performed for 3D slope stability(e.g.Hicks and Spencer,2010;Hicks and Li,2018)as well as square footing(Ahmed and Soubra, 2012; Al-Bittar and Soubra, 2014; Kawa and Pu?a,2020). To the best of the authors’ knowledge, in the case of strip foundation in three dimensions, there are only two such studies(Simoes et al., 2014; Kawa and Pu?a, 2020) which are based on random finite element limit analysis (RFELA) and random finite difference method(RFDM),respectively.In both of these works,the high computation cost of such analysis is highlighted. This cost is much higher when it is necessary to estimate small probabilities(since a large number of simulations need to be performed). This situation indicates the need of finding more efficient approaches that account for 3D spatial variability but at the same time allow to perform probabilistic analyses in a reasonable amount of time.

An efficient algorithm for accessing the upper bound of bearing capacity in the case of two-layered soil based on RFMM mentioned above, which (at least partially) accounts for the 3D spatial variability, was shown in Chwa?a and Pu?a (2020). A plain failure mechanism proposed by Micha?owski and Shi (1995) was used to analyse the bearing capacity in the case of the following arrangement of layers in the subsoil: a relatively thin, strong, frictional layer lying on top of a soft cohesive layer. This arrangement corresponds to the typical problem of a working platform (made of sand) which lies on a layer of soft clay (BR 470, 2004). Since the upper layer corresponded to a well-controlled man-made structure, the values of parameters were assumed as constant and the spatial variability concerned only the lower cohesive layer modelling natural soil. Despite the plain failure mechanism used in this layer, 3D spatial averaging was performed. This joined approach(2D mechanism and 3D averaging) allowed to strongly reduce the computation time needed to assess the probability distribution of foundation bearing capacity.

Another simplification often used when assessing random bearing capacity is the assumption of stationarity of the random field (Doob,1953). A direct consequence of this assumption is the inability of the mean value of the considered parameter to change with depth (the mean values of the parameters have no trend).While for many parameters, this assumption may be correct (c.f.Kawa et al., 2019) in the case of undrained shear strength, the increasing trend in its value with depth is often reported in the literature (e.g. Jaksa and Kaggwa,1996). For this parameter, to the stationary random field, a trend should be added. Usually, a linear trend with the depth is chosen for this purpose(c.f.Li et al.,2015).This issue was omitted in the authors’ previous work(Chwa?a and Pu?a, 2020) where for modelling of undrained shear strength, a stationary field was used. Moreover, as has been mentioned, the proposed method was based on the kinematic approach and represented an upper bound assessment of bearing capacity.In such a situation, the question arises as to how far this upper assessment differs from the exact solution. To answer this question, it is necessary to verify the results using another method for estimation of the exact solution. In the case of probabilistic models, such verification must concern not only the mean values but also the variance of the solution. The discussed study by Chwa?a and Pu?a(2020)did not include such verification.

The present work is an extension of the work by Chwa?a and Pu?a (2020). As before, for the upper bound assessment of the bearing capacity of the working platform,a joined approach is used,i.e. for the assumed plain failure mechanism, the 3D spatial variability is accounted for by averaging of the undrained shear strength (which is the only random parameter) performed also in the out-of-plane direction. However, the description of undrained shear strength for the natural clay layer was modified by adding a linear trend with depth to the stationary random field. Different values of the slope of this trend, horizontal fluctuation scale, and platform thickness are tested. These results allow for general conclusions about the safety of this type of construction on cohesive soils with undrained shear strength linearly increasing with depth.In addition,for selected cases,the results of the probabilistic upper bearing capacity assessment obtained with RFMM are verified by numerical analysis carried out using RFDM. Two types of analyses are performed for this purpose.Since RFMM is limit state approach,first the associated flow rule is assumed.In this case,the agreement between two approaches is very good, despite the fact that the assumed failure mechanism for RFMM cannot change its geometry along the foundation (the 3D analysis only accounts for 3D variability)while this is possible using RFDM approach.Moreover, the efficient RFMM method always provides a conservative estimate of RFDM results. For this reason, in case of associated flow rule, it seems to be preferable method. In the second case, a nonassociated flow rule is assumed for sand layer. As shown in this case, the results based on the RFMM provide a non-conservative estimate for the mean value and a conservative estimate of the coefficient of variation for the values obtained using RFMM.Based on those results,for non-associated flow rule,a combined approach utilizing the results of both methods to provide a conservative estimate of structure failure is proposed.

2. Methods

2.1. RFMM: failure mechanism and spatial averaging

In the present paper, the problem of bearing capacity of the working platform is analysed using RFMM developed and later modified by one of the authors (Pu?a and Chwa?a, 2015, 2018;Chwa?a, 2019, 2020; Chwa?a and Pu?a, 2020). This approach employs the Monte Carlo framework.To assess the bearing capacity of strip foundation in individual Monte Carlo realisation, a kinematic method utilizing upper bound limit state theorem is used. The failure mechanism assumed in the study for two-layered subsoil was adapted from study by Micha?owski and Shi(1995).However,a special version of this mechanism,accounting for spatial variability of soil strength parameters(modelled by a random field)proposed by Chwa?a and Pu?a (2020), is used in the numerical analyses. The considered failure geometry is shown in Fig.1.This mechanism(as in the previous work of the authors) is used to model the arrangement of soil layers,which is typical for the case of a working platform.The upper layer,made of sand,as an anthropogenic layer formed in well-controlled conditions is assumed to be made of homogenous (not variable), frictional Mohr-Coulomb material;thus,the velocity jump vectors in this layer are inclined at the same angle equal to internal friction angle φ (see vectors v1-v5, v10-v15and v20in Fig.1).The bottom layer of clay is assumed to be made of natural (spatially variable) and undrained cohesive (following Tresca criterion) material. Due to the spatial variability of undrained shear strength(this parameter is assumed as the only one modelled by random field) considered in this layer, the failure mechanism is non-symmetrical.As shown in Fig.1,the mechanism used consists of 11 rigid blocks.Two blocks numbered 6 and 11 are located only in the sand;the others consist of two soil layers.Note that (as is the original mechanism by Micha?owski and Shi,1995)there is no velocity discontinuity on the line separating the two layers. Additionally to the exemplary geometry of the failure mechanism, the velocity hodograph is also shown in Fig.1. Please note that in the figure, 3D failure mechanism is presented. This is the mechanism considered here for spatial averaging;the length of the domain(length of averaging)in the y direction is denoted as a.However, since the considered mechanism is only an extrusion of plane mechanism (its geometry along the foundation is constant),the formula used for bearing capacity is for plane-strain condition,the dimension a influences only the area of averaging of undrained shear strength in the bottom layer of subsoil.According to Chwa?a and Pu?a (2020), the bearing capacity can be calculated as where Fjdenotes the area of rigid block j(see Fig.1),liis the length of dissipation surface i, γ is the unit weight of soil, and viis the magnitude of the velocity jump vector along the dissipation line i.Note that[vi]denotes the vertical component of velocity vi(negative if upward). The undrained shear strengths that correspond to dissipation lines liare denoted by ci.More details on the derivation of Eq. (1) can be found in the study by Chwa?a and Pu?a (2020).According to the cited work, the assumed failure mechanism consisting of 11 blocks provides a compromise between accuracy and numerical efficiency.

Eq. (1) provides the value of bearing capacity for the assumed failure mechanism geometry. However, to find the lowest value of bearing capacity for the assumed mechanism, its geometry needs to be optimized; to find this optimal geometry, a proper optimization procedure needs to be applied. In this study, a procedure based on the simulated annealing (Kirkpatrick et al.,1983) is used and following parameters are subjected to modifications within the optimization procedure: xP1，…，xP10，xC1，xC2，l2，l3，l4，l10，l11，l12(details are provided in the numerical algorithm section).

To use Eq. (1), the spatially variable undrained shear strength needs to be averaged over dissipation regions in the cohesive layer.The algorithm of spatial averaging is based on Vanmarcke(1977a,b,2010)’s spatial averaging concept.This approach does not require a direct generation of a random field X. Instead of that, a vector of correlated random variables is generated in such a manner that each component of vector XVcorresponds to the dissipation region resulting from the failure mechanism geometry. The general formula for the spatial average is given as

where X is the initial random field that describes undrained shear strength spatial variability and Viis the averaging domain (the considered dissipation region). Note that due to the stationarity assumption(the trend will be added later),only variances of the XVcomponents are being reduced(the mean values are preserved),i.e.Var(XV，i) = γi(V)σ2X, where γi(V) denotes the variance reduction function. As shown in Fig. 2, the dissipation regions considered in this study are rectangles.The exemplary formulae for variance and covariance between two rectangles are given at the end of thissection. The auto-correlation function for spatially variable layer was assumed to be of separable Markovian type:

Fig. 2. Dissipation regions used for covariance matrix determination. For the two distinguished rectangles, the variances and covariance are derived in the text.

where θx,θyand θzare the fluctuation scales (Vanmarcke,1977a;Fenton and Griffiths, 2008). Despite distinguishing fluctuation scales in three directions, in the following sections, values of both scales for horizontal directions are identical and denoted as θh( = θx= θy). By analogy, θzis denoted as vertical fluctuation scale θv. The covariance function can be obtained from Eq. (3) by multiplying the right-hand side of Eq. (3) by the variance of the initial random field σ2X, i.e. R(Δx，Δy，Δz) = σ2Xρ(Δx，Δy，Δz). For 11 rigid blocks assumed in the mechanism(Figs.1 and 2),there are 16 dissipation regions for which the new variances, as well as covariances between each pair, should be determined. Therefore,the resulting size of the covariance matrix is 16 ×16. Calculation of the covariance matrix components requires dedicated formulae. Their derivation is possible through the usage of the following equation (see Pu?a and Chwa?a, 2015):

By using Eq. (4), the covariance between random variables XViand XVjcan be determined. As an example, formulae for variances and covariance for the two averaged dissipation regions (two rectangles) associated with block 2 distinguished in Fig. 2 are derived below.

The variance derivation for the dissipation region P1A1A′1P′1needs a parametric representation of this region. The following parametrization is introduced:

By substituting Eqs. (3) and (5) into Eq. (4), the following formula for the new variance (after averaging) is derived:

The covariance corresponding to the two above discussed dissipation regions can be calculated by

Other variances and covariances are derived by analogy to those shown above. Finally, a covariance matrix can be determined (Eq.(10)).The indices in the covariance matrix(as well as in Eqs.(6),(8)and (9)) correspond to the indices of length lishown in Fig.1.

The resulting covariance matrix CXis symmetric and positively definite.Due to the symmetry of the covariance matrix,only 136 of the total 256 components need to be computed.The derived matrix CXis used in the procedure of averaged and correlated undrained shear strengths generation that is described in the next section.

2.2. Numerical algorithm of RFMM

In Chwa?a and Pu?a (2020), a stationary random field was considered.However,in the present study,the algorithm proposed by Chwa?a and Pu?a (2020) has been extended to allow for the consideration of the linear trend in undrained shear strength in the spatially variable soil layer. Additionally, instead of a Gaussian correlation function,a Markovian one is used(see Eq.(3)).A newly proposed algorithm is described below. For a description of the spatially variable clay layer, a lognormal random field is used.Lognormal random fields are often reported in the literature to be used for soil strength spatial variability description(e.g.Fenton and Griffiths,2008).The values following this distribution can be easily obtained by the exponential transformation of values following the so-called underlying normal distribution. The same procedure is used in the algorithm below. The correlation structure is assumed to be anisotropic - vertical (θv) and horizontal (θh) fluctuation scales are distinguished.The values of analysed pairs of fluctuation scales are given in the next section;generally,greater values for θhthan for θvare assumed. As have been mentioned, the simulated annealing is used to optimize the geometry of the failure mechanism.Please note that in the below algorithm,simulated annealing procedure is used twice for one Monte Carlo realisation. The first simulated annealing call is focused on finding the optimal geometry assuming that the variables associated with individual slip lines are independent.The result of this first call (the first optimal geometry found) is used to calculate the covariance matrix. Then the simulated annealing is called for the second time.This time the optimal geometry is found using the covariance matrix obtained before. In this way, the values of variables in the second step of simulated annealing are properly averaged (their variance is reduced) and correlated. The subsequent steps of the proposed algorithm are described in detail in the following paragraphs. The authors have adopted the convention of giving details of assumptions inside the description of a given step. Therefore, the descriptions of the individual steps are sometimes extensive.However, to make the reading of the algorithm more clear, a flowchart with general comments has been added and shown in Fig. 3. Both the detailed description of the algorithm and the flowchart are focused on one Monte Carlo realisation;as explained below,all the steps of the procedure(starting from Step 3)have to be repeated an assumed number of times to obtain the bearing capacity distribution:

Fig. 3. Flowchart of the numerical procedure.

(1) Step 1 introduces the values of the initial geometry parameters, i.e. foundation width b, sand layer thickness t, averaging length a (see Figs. 1 and 2), and soil parameters, i.e.friction angle for top layer φ,the mean value μcuand standard deviation σcuof undrained shear strength for the bottom layer,the trend slope αTand fluctuation scales θvand θh.Set i = 1 and start the Monte Carlo loop.

(2) Step 2 continues the Monte Carlo loop: if i ≤ N, follow the steps below; if i>N， go to Step 14 (N is a number of Monte Carlo realisations).

(3) Step 3 generates a random vector Ycuvia a pseudo-random number generator. The generated values are independent and normally distributed with a mean value μY，cuand standard deviation σY，cucorresponding to the underlying normal distribution for the stationary lognormal random field (μcu,σcu) (see Step A in Appendix A). Note that in this step, no additional actions are taken into account to consider the trend in undrained shear strength.For the assumed number of rigid blocks,16 values are generated,i.e.Ycu= (Ycu，1， …，Ycu，16).

(4) Step 4 assumes initial failure geometry vector G =

(xP1，…，xP10，xC1，xC2，l2，l3，l4，l10，l11，l12) (chosen arbitrarily).

(5) Step 5 calculates lognormally distributed undrained shear strengths (by exponentially transforming values of underlying normal distribution) Xcu= (exp(Ycu，1)，…，exp(Ycu，16))?Additionally,in this step,extra part resulting from the linear trend of undrained shear strength is added to Xcu.The vector of these additional values Δcu= (Δcu，1，…，Δcu，16) is calculated based on the locations of the slip lines centres as shown in Fig. 4. For example, according to Fig. 4, Δcu，1=aT(-zA，1-t)/2 and Δcu，5= αT[( -zA，1-zA，2)/2 -t],where αTdenotes the slope of the trend.(The above calculations are performed using the values of coordinates on the global axis z, directed upward from the surface. Alternatively, the local axis zT, starting from the top of the clay layer and directed downward,can be used,which simplifies calculations).Next,a final vector of undrained shear strengths is calculated:cu= Xcu+Δcu.Note that at this stage,the coordinates of Xcuare independent.

(6) Step 6 calculates bearing capacity based on the assumed failure geometry and undrained shear strengths cucalculated in Step 5.

(7) Step 7 initiates simulated annealing based optimization

scheme directed toward minimization of bearing capacity upper bound (first call). The optimization procedure is constructed to take into account the changes in undrained shear strengths resulting from the changes in locations of the centres of the dissipation lines caused by slight variations in the failure mechanism geometry parameters that are introduced within the optimization procedure. It means that during optimization, a vector Xcuremains unchanged(Chwa?a and Pu?a, 2020); however, Δcuare simultaneously updated.This ensures a trend consideration.As a result of the optimization procedure,a failure geometry that delivers the lowest possible bearing capacity upper bound value for the assumed mechanism is found. In Steps 7.1-7.8, a detailed algorithm for subset simulation is given:

(i) Step 7.1 assumes initial failure geometry G from Step 4 and corresponding bearing capacity from Step 6 as the current bearing capacity pc.

(ii) Step 7.2 assumes the control parameters for the optimization scheme of simulated annealing (α， T， Tmin， h),where T is a parameter that modifies acceptance probability (see Step 7.7); Tminis the minimum value of T at which the simulation ends(see Step 7.3);α is responsible for reducing T, i.e. T = αT as shown in Step 7.8; and h is the number of simulation for each value of T.According to the earlier experiences(Chwa?a and Pu?a,2020)and tests performed for the present case (with trend consideration), the following values are assumed in the study:h = 2000, α = 0?9, T = 10 (starting value), Tmin=10-8.

Table 1 The values of parameters used for RFMM.

(iii) Step 7.3 checks if T >Tmin, if true, go to Step 7.4 and set j = 1; if not true, end algorithm and assign as optimal failure geometry a vector Gopt= G and popt= pcas corresponding optimal bearing capacity. Next, go to Step 8.

(iv) Step 7.4 checks if j ≤h,if true,go to Step 7.5;if not true,go to Step 7.8.

(v) Step 7.5 creates a so-called neighbouring set of failure geometry parameters to G by introducing slight random modifications of G. A new vector Gnis obtained by the following modification: Gn= G+ ΔG, where ΔG =(U1[- 0?01 m， 0?01 m]， …， U18[- 0?01 m， 0?01 m]), in which U[·] is a uniform distribution on a given interval.

(vi) Step 7.6 calculates the bearing capacity for the neighbouring set of failure geometry Gn, and the resulting bearing capacity is denoted by pn(candidate bearing capacity).

(vii) Step 7.7 updates the value of current bearing capacity,i.e.if pc≥pn, update pcas pc= pn; if pc

(viii) Step 7.8 updates T = αT and go back to Step 7.3.

(14) Step 14 ends the numerical procedure.

As mentioned above, the value of bearing capacity obtained in Step 13 of the presented algorithm is assumed as the final for the considered Monte Carlo realisation.Please note that the correlation matrix used for this result (obtained in Step 8) corresponds to the geometry of failure mechanism obtained in Steps 7.1-7.8 which differs from the final geometry of the failure mechanism obtained in Step 13. However, as shown in Chwa?a (2018), subsequent updating of the correlation matrix (using the final geometry) has negligible influence on the result.

The capability of the presented algorithm to account for the linear trend in cuwas first tested for non-random cases. All the values of parameters used are summarized in Table 1. The results showing the influence of trend slope on the size of the mechanism for the considered cases is presented in Fig.5.The figure shows the optimal failure geometry for the case of no trend in the bottom cohesive layer(αT= 0)and for three values of the trend slope,i.e.αT= 9 kPa/m, 18 kPa/m and 36 kPa/m (the assumed sand layer width is 0.5 m).It is worth mentioning that the greatest differences in the failure geometry size are observed between cases of αT= 0 and 9 kPa/m. As shown in Fig. 5, trend consideration provides a significant reduction in the failure mechanism size. Intuitively, in the case of a lower trend slope, the failure zone is deeper and the slip lines are longer.Hence,the averaged values of cuon individual lines differ from each other(due to different z coordinates of their centres).In the case of a larger slope of the trend,the failure zone is shallower and the averaged values of cuare more uniform than in the previous case. The details on this issue for the case of probabilistic modelling are discussed in Section 4.

Fig. 4. Illustration of linear trend consideration in bearing capacity calculations.

Fig.5. Failure geometry in the deterministic case(undrained shear strength modelled with constant value of 9 kPa with trend) for a top sand layer width of 0.5 m. The obtained failure geometries are symmetrical.

Fig. 6. Exemplary failure geometry obtained for spatially variable undrained shear strength in the bottom layer (αT = 9 kPa/m, θv = 1 m and θh = 2 m).

Scenarios shown in Fig. 5 are obtained for the deterministic case; thus, the failure geometries are symmetrical respect to the foundation centre. In a random case, if the spatial variability of undrained shear strength is considered, the failure geometries are nonsymmetric (since the soil on one side of the foundation is always weaker).An example of such mechanism obtained for a single realisation of the field of undrained shear strength, for αT=9 kPa/m, θv= 1 m and θh= 2 m, is shown in Fig. 6. Except the field being random, the same values of parameters as in deterministic case were used(Table 1).Note that the results presented in the figure only concern one realisation. To obtain the mean values and coefficients of variation for the bearing capacity, at least hundreds of realisations are needed(see number N in Step 2 in Fig.3).As mentioned above, the method is quite efficient, which is preserved when adding a linear trend.The average computation time for a 4-core i7 processor takes about 20 s per one realisation(about 1.5 h for 300 realisations). What is important, this value is independent of the assumed averaging length a.

It is worth recalling that the method used is based on the kinematic method, and by the limit state theorem, represents an upper estimate.To verify the quality of this estimate,apart from the analyses incorporating this approach,calculations using the RFDM were also carried out in this study.

2.3. RFDM based verification

The RFDM approach used to verify the results obtained using RFMM also employs the Monte Carlo framework. To assess the bearing capacity of a strip foundation on the working platform,an individual realisation of the boundary value problem was modelled in FLAC3D, a finite difference method code. Few series of computations with different values of the modelled domain length a(the length of the modelled domain corresponds to the averaging length in RFMM), trend slope αT, and horizontal fluctuation scale θhwere performed.Regardless of the length of the domain,a foundation in the direction parallel to its axis was modelled along the entire length of the domain. In a plane perpendicular to its axis, the foundation was located at the centre of the domain (due to the random nature of the field, the problem, as in RFMM, was not symmetrical). The boundary conditions were assumed as follows:for each side surface of the domain,either perpendicular or parallel to the foundation axis, displacements in the direction perpendicular to the surface were fixed, and on the bottom surface, all displacement directions were fixed.

In the simulation, a foundation width b of 1 m was assumed(identically as in the reference problem using RFMM). The thickness t of the top friction layer modelling the working platform was assumed to be 0.5 m. An elasto-perfectly plastic model of the material was assumed for both soil layers.The bearing capacity of the foundation was determined by controlling the vertical displacement of the nodes on the surface of the platform (modelling the foundation), determining the sum of the reactions at these nodes,and then dividing it by the foundation area. This value gradually increased as the nodes were pressed against the soil until finally a failure mechanism was established for which a stable value was obtained. No additional ties were imposed on the horizontal movements of these nodes (smooth foundation). A node vertical velocity (directed downward) of 5 ×10-6m per computation step was assumed,and a stable solution was obtained after 10,000 steps.These values were determined by trial and error testing that for all analysed cases (random field realisations), the final result was always obtained without loss of numerical stability.

For the top layer of sand (which, as with the RFMM, was assumed to be non-random), the model parameters are listed as follows: Young’s modulus E1= 50 MPa, Poisson’s ratio υ1= 0.25,friction angle φ1=35°and unit weight γ=18 kN/m3.The values of parameters correspond to RFMM (if the parameters exist in both methods) and are typical for medium to dense sand materials;comparable values of parameters were reported as measured for the working platform material using DMT method in Balachowski and Bialek (2015). It should be noted that the values of elastic parameters should have negligible influence on the results (c.f. Pu?a and Zaskórski, 2015). Since the RFMM is a limit state method, an associated flow rule was assumed resulting in the value of dilation angle ψ1=35°.The influence of the flow rule was later checked by repeating the simulation for ψ1=0°.The 0°value is lower than that obtained using the classical formula ψ=φ-30°,but the difference is not very large and it was chosen as the most conservative value.In all simulations, for computation reasons, a small value of cohesion c1=1 kPa for the platform material was also assumed (the same value of parameters including cohesion was assumed for the comparative analysis performed using RFMM). In the case of the clay layer,the elasticity parameters and unit weight were assumed as non-random with values E2=20 MPa,υ2=0.48 and γ=18 kN/m3. In addition, two plasticity parameters, i.e. the friction and dilation angles, were assumed as non-random and 0°. Undrained shear strength was modelled by a stationary random field (with correlation function expressed by Eq. (3)) with the added linear trend increasing with the depth. Mean values for the stationary field was assumed as μcu= 9 kPa and standard deviation as σcu= 4.5 kPa(coefficient of variation equals 0.5). All these parameters are summarized in Table 2.

Individual realisations of the stationary part of the field were generated using a 3D version of the Fourier series method (FSM)(Jha and Ching,2013b).The vertical scale of fluctuation θv=θzwas assumed as 0.5 m. Two values of horizontal scales of fluctuation were considered, i.e.θh=θx=θyequal to either 2 m or 10 m. The values of the field were assumed to follow a lognormal distribution(identically as in the case of RFMM). However, the local averages generated with FSM followed the underlying normal distribution.The generated values were subsequently transformed exponentially to obtain a lognormal random field (the formulae for the values of parameters of such distribution,based on the parameters of lognormal distribution, were also used in RFMM and are included in Step A in Appendix A). The small influence of such transformation on the values of fluctuation scales (Huang and Griffiths, 2015; Pu?a and Griffiths, 2021) was omitted. Adding the trend values dependent on the individual finite difference zone centroid was the last step of the field generation. The considered values of trend slope αTwere 9 kPa/m and 18 kPa/m.

Table 2 The values of parameters used for RFDM.

Fig.7. Geometry of the model and values of the parameters:(a)Model dimensions and deterministic value of friction angle, and (b) Exemplary realisation of the field of undrained shear strength including the trend for θx = θy = 2 m (θz = 0.5 m) and trend slope αT = 18 kPa/m.

The four problems mentioned above (two values of the horizontal fluctuation scale θh×two values of the slope trend αT)were modelled for different values of length of the domain, i.e. 0.167 m(one zone),1 m, 2 m, 4 m and 8 m. The geometry of the problem together with exemplary realisation of the undrained shear strength for trend slope αT=18 kPa/m and θh=2 m(θv=0.5 m)is shown in Fig. 7. The nodes associated with the foundation are marked with vertical black lines. As seen, the cohesion (in the second layer corresponding to the undrained shear strength),which in the first layer equals 1 kPa, in the second layer is increasing linearly but with random residuals.

Domain size in the direction perpendicular to the foundation axis was adopted to avoid the plastic zone reaching the domain boundaries (cf. Kawa and Pu?a, 2020). The necessary domain size satisfying this condition for all random cases considered was determined as 25.5 m for width and 3.5 m for height(including the sand layer). For cubic zones with all dimensions equal to 0.167 m used to discretize the domain,this corresponds to 153×21 zones in one slice (0.167 m along the foundation axis). The relatively large width of the domain considered is due to the wide range of plastic zones in the sand layer (when modelling a smaller domain width,the plastic zones at the surface reached the edge of the domain).The large volume of domain resulted in a rather long computation time. For example, in the case of a modelled foundation length of 8 m, the computation time for a series of 300 simulations (single value of the horizontal fluctuation scale θh,and a single value of the trend slope αT)for a modern PC processor running on three threads was equal to about two weeks. For that reason, the number of simulations was limited to N = 300. The convergence of Monte Carlo results indicates that this assumption has a limited impact on the obtained accuracy (see Appendix B).

Fig.8. Plastic zones under foundation in exemplary realisations with different values of modelled length a,trend slope αT=18 kPa/m and θx=θy=2 m(θz=0.5 m):(a)a=0.166 m,(b) a = 1 m, (c) a = 2 m, and (d) a = 4 m.

Fig.9. Typical plastic zones under foundation obtained for modelled length a=4 m:(a)αT=18 kPa/m and θh=2 m,(b)αT=9 kPa/m and θh=2 m,(c)αT=18 kPa/m and θh=10 m,and (d) αT = 9 kPa/m and θh = 10 m.

Fig.10. Alternation of depth of the plastic zones along axis of the foundation obtained for a=8 m and θv=0.5:(a)αT=18 kPa/m and θh=2 m,(b)αT=9 kPa/m and θh=2 m,(c)αT = 18 kPa/m and θh = 10 m, and (d) aT = 9 kPa/m and θh = 10 m.

Fig. 8 shows the plastic zones obtained in exemplary realisations for ψ = 0°and different modelled foundation lengths (for ψ=35°,the obtained plastic zones are slightly wider).Note that in the images only the middle part of the model is zoomed(～15 m).The zones which are currently in plastic state are marked with brown colour and the zones which are in elastic state and never passed through the plastic state during the whole process of solving the problem are marked with blue colour.Due to the way FLAC3D handles the problem(iterating while moving the nodes to model the foundation) in individual cells, temporary or constant plastic unloading occurs.This means that the zones,which in the previous numerical time step were in the plastic state,in the next step can be once again in the elastic state.These zones are marked with pink colour as “formerly plastic” (currently elastic) to distinguish them from “purely elastic” blue zones. The results presented were obtained for a fluctuation scale θh= 2 m(θv= 0.5 m) and a trend slope αT= 18 kPa/m. For the modelled foundation length of 8 m (not shown in this figure), the plastic zones obtained were similar in shape and corresponded to those shown in the figure (being the most similar to those obtained for the modelled length of 4 m). No particular influence of the modelled length of foundation on the shape of plastic zones was observed(except that for longer modelled length,the shape could differ along the foundation).

The use of a full 3D simulation allows the failure mechanism to be alternated along the foundation axis,which is impossible for the assumed failure mechanism in the case of the RFMM.Figs.9 and 10 compare the most typical shapes of plastic zones obtained both in the transverse direction and longitudinally along the foundation for all four cases(two θhand two αT)analysed.Images of the obtained plastic zones showing their shape in the transverse direction are shown in Fig. 9. To demonstrate the alternation of plastic zones depth along the foundation,Fig.10 shows half of the domain cut by a vertical plane parallel to the foundation axis.Half of the domain is presented only to better illustrate the changes in plastic zones along the foundation axis (due to the lack of symmetry of the problem, the numerical analysis involved the whole domain). As can be seen in the figures,for the cases with lower αTvalues(Cases b and d in Figs. 9 and 10, respectively), the plastic zones (corresponding to the failure mechanisms)reach deeper.Moreover,for a smaller value of the horizontal fluctuation scale θh= 2 m (Cases a and b in Figs. 9 and 10, respectively), plastic zones, which are asymmetrical and more frequently change with the length of foundation than in the case of the considered longer horizontal fluctuation scale, were obtained.

3. Results of numerical analyses of RFMM

Examination of the impacts of the trend slope αT, averaging length a, homogenous sand layer thickness t, and vertical and horizontal fluctuation scales on the probabilistic characteristics of bearing capacity is performed in this section. For the analyses, the RFMM procedure described in the previous sections was used.Both of mean value and coefficient of variation for bearing capacity were estimated based on 1000 Monte Carlo realisations.In all considered cases, it was assumed that the stationary random field was described by the mean value of undrained shear strength μcu=9 kPa and coefficient of variation vcu= 0?5.As detailed earlier,for the stationary random field, an additional part was added that results from the trend assumption. Due to the assumed undrained conditions in the spatially variable layer,no correlation with other soil parameters is assumed in the study. Moreover, the results obtained in this section were obtained under the assumption of small cohesion in the homogenous sand layer, i.e. c = 1 kPa. First,the impact of the trend slope of cu,i.e.a value of αT,was examined.The results obtained for θv= 1 m and 0?5 m are shown in Figs.11 and 12, respectively. As shown in Fig. 11a and b and 12a and b,the mean value μ of bearing capacity increases with an increase in trend slope αT. The situation is not so clear for the coefficients of variation v of bearing capacity shown in Fig.11c and d and 12c and d,which exhibit interesting behaviour especially visible in Fig.12c and d, i.e. a local maximum is observed. This can be explained by two opposite phenomena. On the one hand, the coefficients of variation are expected to decrease with increasing trend slope αT.This is because the greater part in cudepends on the trend(which is non-random) and not on the random field residuals. Moreover, a greater slope of the trend results in a shallower failure geometry as shown in Fig. 5, and the resulting percentage of bearing capacity that originates from the homogenous sand layer slightly increases.On the other hand, the shallower failure geometry results in smaller slip lines in the cohesive layer, which in turn results in smaller averaging areas.Because of that,the resulting variability of a random part of cu，avgbecomes greater.Although there seem to be more effects in favour of decrease of the coefficient of variation with increasing trend slope (greater influence of the homogenous sand layer and larger non-random part of cu，avg),which is reflected by most of the results(Figs.11 and 12,all graphs for αTgreater than 9), a small local increase in the coefficient of variation resulting in the mentioned local maximum is also possible.

Fig.11. Mean values and coefficients of variation for bearing capacity as functions of the trend slope αT:(a,c)a=1 m,and(b,d)a=8 m.Vertical fluctuation scale is assumed in all cases as θv = 1 m.

Fig.12. Mean values and coefficients of variation for bearing capacity as functions of the trend slope αT:(a,c)a=1 m,and(b,d)a=8 m.Vertical fluctuation scale is assumed in all cases as θv = 0.5 m.

Fig.13. Mean values(a)and coefficients of variation(b)for bearing capacity as functions of horizontal fluctuation scale θh for different values of trend slope αT,vertical fluctuation scale θv and foundation length a.

The dependencies of μ and v on horizontal fluctuation scale θhfor two vertical fluctuation scales θv= 0?5 m and 1 m, two trend slopes αT= 9 kPa/m and 18 kPa/m and two averaging lengths a =1 m and 8 m are shown in Fig.13.As shown,the trend slope has the greatest impact on bearing capacity mean values among the considered cases. For each case, higher μ are obtained for greater vertical fluctuation scale (θv= 1 m) and for longer foundations(a = 8 m). Bearing capacity mean values generally decrease with an increase in θh;however, for over half of the considered cases, a local minimum near θh= 10 m is observed.This can be interpreted as a weak worst-case effect occurrence in mean value (e.g. Ching et al., 2017; Pieczy′nska-Koz?owska et al., 2017). As expected, for increasing θh,greater values of v are observed in all cases.However,the rate of increase for relatively large θhis greater for foundation length a = 8 m in comparison with a = 1 m. This is because for a = 1 m,the size of the failure mechanism is relatively small when comparing with θh(further increase in θhdo not change significantly the resulting variability of cu). The effect of foundation length influence on the observed v is even better visible when comparing relative differences in v for two considered foundation lengths (a = 1 m and 8 m) for the smallest θh= 1 m with the greatest θh= 20 m.In the case of θh= 1 m,those differences are about 100%; but, for θh= 20 m, those differences are reduced to approximately 10%-15%. This example indicates the importance of considering 3D spatial variability and the actual value of foundation length for analysis of strip foundation when spatial variability of the soil is considered. It also shows the importance of precise determination of horizontal fluctuation scale which in many cases is problematic (e.g. Cami et al., 2020).

The dependence of the mean values and coefficients of variation for bearing capacity on the thickness of the homogenous sand layer(t)and horizontal fluctuation scale is shown in Fig.14.As expected,the mean value of bearing capacity increases if a case with greater t is considered. The inverse is observed for v, where for greater t,lower values of v are observed. This effect illustrates the ability of the homogenous sand layer to suppress the impact of the soil spatial variability in the bottom layer.

Table 3 Scenarios considered in the analyses performed by RFMM.

The above results indicate several important issues. The first is to take into account the 3D soil spatial variability even if we consider long foundations for which, in a deterministic case, a plane-strain condition is commonly assumed. If the foundation is rigid (e.g. it is stiffened by a relatively high wall), its length significantly influences the coefficient of variation of bearing capacity:its value decreases with the increase in foundation length.The second issue is the trend consideration in cohesive soil. As shown, the trend significantly affects the bearing capacity probabilistic characteristics and ignoring it by assuming a stationary random field may result in strongly inaccurate values, which is even more important in the case of probabilistic modelling. The intuitive behaviour of the coefficient of variation for bearing capacity is shown in Fig.14,where the impact of natural soil spatial variability decreases with the increase in homogenous layer thickness. To summarise the numerical analyses performed by RFMM, all scenarios considered in this section are listed in Table 3.

Fig.14. Mean values (a) and coefficients of variation (b) for bearing capacity as functions of horizontal fluctuation scale θh for different widths of homogenous sand layer t under αT = 9 kPa/m, θv = 1 m and a = 1 m.

4. Comparison between RFMM and RFDM

Comparison of the results of the approach based on the RFMM and RFDM was performed for the cases analysed using RFDM, i.e.two trend slopes (αT= 9 kPa/m and 18 kPa/m), two pairs of fluctuation scales (θv= 0.5 m and θh= 2 m, and θv= 0.5 m and θh= 10 m), and only one value of the thickness of the sand layer(t=0.5 m).The assumed values of slope in the trend are taken in a frequently reported scope. In addition, the assumed values of fluctuation scales were taken from the range reported in the literature to cover both the cases of relatively small anisotropy and large anisotropy. The RFMM results were obtained for N = 1000 Monte Carlo realisations.Due to the significant computation time required for RFDM, the number of Monte Carlo realisations for each considered case was reduced to 300 (for the same reason, only selected cases were analysed with this method).As described in the numerical algorithm section, a 3D soil spatial variability is considered in both methods.For this reason,a parameter a(in the case of RFMM, averaging distance and in the case of RFDM, modelled length of foundation) is necessary to be assumed. The following values of a are investigated in the comparison of both methods:a = 0?167 m, 1 m, 2 m, 4 m and 8 m. The a = 0?167 m corresponds to the size of the assumed mesh in RFDM and thus it is the smallest possible value for which the methods can be compared.

4.1. Associated flow rule

Since the RFMM is a limit state method, the assumption of the associated flow rule for RFDM(ψ1=35°for sand)seems to be the most appropriate for comparison. The results obtained for both methods in terms of mean value μ of bearing capacity and the corresponding coefficient of variation v are shown in Figs.15 and 16 for αT= 9 kPa/m and 18 kPa/m, respectively.

As shown in Figs.15 and 16,the mean values of bearing capacity for both approaches are in good agreement for all considered scenarios.This agreement can be surprising taking into account the mentioned difference between the failure shapes considered in both methods; only in the RFDM, the failure shape along foundation axis has been changed. Despite this fact, the mean values obtained using RFDM for two different fluctuation scales are actually almost identical. The small differences visible in particular cases can be the effect of averaging performed on relatively small number of realisations. As it appears from the convergence analysis presented in Appendix B, this is the case for the mean value obtained for αT=9 kPa/m,θh=2 and a=2 m(where the difference from the mean value obtained for other considered fluctuation scale seems greater than that for other values of a presented in Fig.16).Since the same realisations of stationary random fields were used for both considered trend slopes, this is probably also the case for αT= 18 kPa/m, θh= 2 and a = 2 m (for the latter case, the convergence analysis was not preformed).Greater and more stable difference occurs between the values obtained for different fluctuation scales with RFMM. As a result, the values obtained with RFMM for θh= 2 m is usually grater compared to RFDM(which is natural since RFMM is upper bound approach), whereas values obtained for θh= 10 m are lower. The differences are, however,very small and can be the effect of numerical errors.

Greater difference between the results obtained with the two methods can be observed for the coefficient of variation.The values are still in good agreement;however,it appears that in all analysed cases, the values obtained using RFMM are greater. From the analysis of convergence presented in Appendix B, it appears that the almost negligible difference in the results obtained with the two approaches for αT= 9 kPa/m, θh= 2 and a = 0.167 m for a greater number of realisations would be greater, while the relatively large difference obtained for αT=9 kPa/m,θh=2 and a=1 m would be smaller.

Fig.15. Mean values (a) and coefficient of variations (b) for bearing capacity obtained using RFMM and RFDM (associated flow rule assumed for both methods) in the case of αT = 9 kPa/m and θv = 0.5 m and two horizontal fluctuation scales θh = 2 m and 10 m.

Fig.16. Mean values (a) and coefficient of variations (b) for bearing capacity obtained using RFMM and RFDM (associated flow rule assumed for both methods) in the case of αT = 18 kPa/m and θv = 0.5 m and two horizontal fluctuation scales θh = 2 m and 10 m.

It is worth noting that if we treat the RFDM results as potentially closer to a close solution(except the mentioned numerical errors),from the point of view of foundation reliability, the mean values obtained from the RFMM should be considered as precisely but not always conservatively estimated (in fact, as mentioned, we would expect the mean values from RFMM to be always grater) and the values of the coefficient of variation as conservatively estimated.However, taking into account that in structural reliability, we usually deal with very low probabilities,especially in the case of the ultimate limit state, and that the differences between the mean values of bearing capacity estimated using the two methods are negligible,and differences between the coefficients of variation are relatively large, it can be concluded that RFMM can be considered as conservative approach of reliability of the considered structure.This was checked by the authors for all scenarios discussed in this section by a lognormal probability distribution approximation to the resulting bearing capacities. The resulting lognormal distributions were used for calculating the probability of failure and the corresponding reliability index β. For all scenarios, RFMM (due to greater coefficient of variation for bearing capacity) is more conservative than RFDM for reliability index β greater than 2. This means that in the case of reliability analysis in civil engineering(where according to ISO 2394:2015 (2015), the indices are greater than 2), the RFMM can be indeed considered as a conservative approach. Since the analysis of bearing capacity using RFMM are additionally much faster, it seems that in the considered case, this method can be a better choice for reliability assessment at least in the case of associated flow rule. Although the design values obtained with the method can be less precise than those obtained using RFDM or RFEM(for appropriate number of realisations),they are also less prone to numerical errors. Good agreement between coefficients of variation estimated with the two approaches on one hand and efficiency and rigorous conservatism of the RFMM method on the other hand are obvious advantages of this approach.However, the concept of using only this method requires further investigation.

4.2. Impact of dilation angle ψ

The associated flow rule is known to be hardly confirmed by laboratory test results and, maybe more importantly, to provide a non-conservative estimation of the critical load.For this reason,in numerical calculations, a non-associated flow rule with much lower dilation angle than friction angle φ is typically employed.To investigate the impact of dilation angle ψ on the comparison of RFDM with RFMM, the cases analysed with RFDM were computed once again assuming for the sand layer ψ1=0°(which is the most conservative case).The comparison of these results with previously presented results of RFMM is shown in Figs.17 and 18.

Fig.17. Mean values(a)and coefficient of variations(b)for bearing capacity obtained using RFMM(associated flow rule)and RFDM(dilation angle ψ=0°)in the case of αT=9 kPa/m and θv = 0.5 m and two horizontal fluctuation scales θh = 2 m and 10 m.

Fig.18. Mean values(a)and coefficient of variation(b)for bearing capacity obtained using RFMM(associated flow rule)and RFDM(dilation angle ψ=0°)in the case of αT=18 kPa/m and θv = 0.5 m and two horizontal fluctuation scales θh = 2 m and 10 m.

As seen,for both considered slopes of the trend,the assumption of a much lower ψ in RFDM results in a significant drop of the mean value of bearing capacity obtained with this method. According to the results,the difference between mean values obtained with the two considered approaches rises to about 10% (in all cases, the values obtained with RFMM are greater). On the other hand, the change of ψ only slightly increases the coefficients of variation obtained using RFDM.For this reason,the values obtained using the two approaches are in even slightly better agreement than those in Section 4.1 for the associated flow rule. Still, the coefficients of variation obtained with RFMM in all considered cases are greater.The only exception is the value for αT= 9 kPa/m, θh= 2 and a=0.17 m.However,based on the convergence analysis presented in Appendix B,it can be concluded that this is just the effect of small number of simulations for RFDM.

It should be noted that in this case(i.e.for ψ=0°),if we assume that the RFDM results are potentially closer to a close solution, the mean values obtained with RFMM should be considered as nonconservatively estimated and the coefficient of variation as conservatively estimated.Due to relatively large differences in mean values obtained with the two approaches for the assumed value of dilation angle for RFDM as 0°,the application of RFMM is likely to result in a non-conservative estimate of the design load value for the reliability indices considered in this case(ISO 2394:2015,2015).Thus,although this approach is effective, if the non-associated flow rule is considered,it should not be directly used.On the other hand,the well and conservatively estimated coefficients of variation by this method seem to encourage a combined approach that utilises the results of both methods. Correct estimation of the mean value requires much less computation effort(smaller number of realisations)than that of the coefficient of variation.Thus,to assess the mean value of bearing capacity,the results of RFDM probabilistic modelling can be used for a reduced number of simulations (i.e.100-200), while to assess the coefficients of variation, RFMM can be used which is efficient and conservative. However, application of this concept requires further analyses and detailed testing in a wide range of soil and geometry parameters. Note that the differences observed between the two methods in terms of the coefficient of variation(for both associated and non-associated flow rules) are small concerning other possible uncertainties we are struggling with (e.g. in situ determination of horizontal fluctuation scales),which seems to be a positive aspect in the context of the possibility of applying both proposed approaches.

5. Conclusions

In the present work, RFMM was used for the stochastic assessment of the bearing capacity of the working platform on cohesive soil. The modified algorithm which considers not only the spatial variability of the undrained shear strength of the cohesive layer but also the linearly increasing trend for that parameter was presented.The algorithm was then used to assess the influences of trend slope,fluctuation scale as well as the working platform thickness on the probabilistic properties of the bearing capacity of the platform.Finally, the RFDM verification of RFMM results was provided. The following conclusions can be drawn from the paper:

(1) Estimating the bearing capacity of a strip foundation on spatially varying soil is an important and actual issue.Although the problem appears to follow plane-strain conditions, in the case of spatial variability of the soil parameters, it is necessary to analyse the structure in 3D space,which is usually associated with high computation cost. For this reason, estimation of the probability of exceeding the ultimate limit state at an acceptable level, using typical numerical methods (RFEM or RFDM), is usually highly inefficient. An alternative approach could be the use of RFMM.This method is efficient and allows to account for the 3D spatial variability of soil parameters.

(2) Compared with the previous work of one of the authors,the modified RFMM presented here allows, in the case of the analysed cohesive soil, to account for both the spatial variability of the undrained shear strength and the linear trend with depth for this parameter. In the case of analysed cohesive soil in undrained condition,this modification has a major impact on the results obtained. As it is shown, the computation cost of the proposed modification has virtually no impact on the efficiency of the method.

(3) The modified RFMM algorithm was used to show the effects of quantities such as trend slope,horizontal fluctuation scale and platform thickness on the mean values and coefficients of variation of bearing capacity. Each time the results were presented for N = 1000 Monte Carlo realisations. As expected, increasing the trend slope results in an increase in mean value of bearing capacity.Less intuitive but explainable issue is the local maximum occurring in a few cases for the coefficient of variation for bearing capacity, shown as the trend slope function. It is also shown that increasing the horizontal fluctuation scale causes a slight decrease in mean value and a significant increase in coefficient of variation of bearing capacity. Finally, it was shown that increasing the thickness of the working platform layer increases the mean value and decreases the coefficient of variation of platform bearing capacity.

(4) Another novelty in the present work is the quantitative verification of the RFMM approach, being an upper bound estimate of the exact solution, carried out using RFDM. Due to the high computation cost, the number of simulations N for RFDM was limited to 300. The verification results obtained with different lengths of the modelled foundation for two typical values of horizontal fluctuation scale and two values of trend slope indicate a fairly good agreement between the probabilistic bearing capacity characteristics obtained from both methods despite the differences between the approaches.It is found that for both associated and nonassociated flow rules assumed for the sand layer, similar coefficients of variation for bearing capacity are obtained by RFDM.However,the mean value of bearing capacity strongly depends on the flow rule assumed. For the associated flow rule,these values are significantly greater than those for the non-associated one.(5) As shown, for the associated flow rule, the mean value of bearing capacity obtained using RFMM is a non-conservative but precise estimate and the coefficient of variation obtained using the RFMM is conservative estimates for results obtained using RFDM.Consequently,the values obtained using RFMM for small values of probability associated with the ultimate limit state provide a conservative estimate of the design value of bearing capacity. In view of the fact that RFMM is a much faster method, it seems that in the case of associated flow rule, it can be a recommended approach for the considered case. On the other hand, for the nonassociated flow rule assumed for the sand layer (ψ = 0°),the mean values obtained by RFMM are non-conservative and not very precise estimates and the coefficients of variation obtained by RFMM are still conservative estimates obtained by RFDM. In this situation, the design values estimated with RFMM are non-conservative and thus the direct use of the method is not recommended. However, to reduce the computation cost in this case,it seems reasonable to use a combined approach that utilises the results of both methods. The results of inefficient RFDM for a reduced number of simulations can be used to determine the mean value of bearing capacity, while the coefficients of variation can be efficiently and conservatively estimated using RFMM.Application of these concepts requires further analyses and detailed testing in a wide range of soil and geometry parameters, which is currently an object of further studies.Declaration of competing interestThe authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

List of symbols

aFoundation length

αTTrend slope in bottom clay layer

bFoundation width

cuUndrained shear strength

EYoung’s modulus

FArea of rigid block

liLength of dissipation surface i

NNumber of Monte Carlo realisations

RCovariance function

RFRandom field

tWidth of sand layer

viVelocity jump vector

Δx,Δy,Δz Distance in the x, y, and z directions

x， y， zGlobal coordinate system

zTLocal vertical coordinate

βReliability index

γUnit weight of soil

θv,θhVertical and horizontal fluctuation scales

θx,θy,θzFluctuation scales in x， y， z directions

μcuMean value of undrained shear strength

υPoisson’s ratio

ρCorrelation function

σcStandard deviation of undrained shear strength

φInternal friction angle

ψDilatancy angle

Var, Cov Variance, covariance A

ppendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jrmge.2021.06.004.