A computational study on the uncertainty quantification of failure of clays with a modified Cam-Clay yield criterion

In this study, the quantitative effect of the random distribution of the soil material properties to the probability density functions of the failure load and displacements is presented. A modified Cam-Clay failure criterion is embedded into a stochastic finite element numerical tool for this purpose. Various assumptions for the random distribution of the compressibility factor κ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa$$\end{document}, of the constitutive relation, the critical state line inclination c of the soil, and the permeability k have been tested and assessed with Latin hypercube sampling followed by Monte Carlo simulation. It is confirmed that both failure load and displacements follow Gaussian normal distribution despite the nonlinearity of the problem. Moreover, as the soil depth increases the mean value of failure load decreases and the failure displacement increases. Consequently, failure mechanism of clays can be determined in this work within an acceptable variability, taking into account the soil depth and nonlinear constitutive relations which in the analytical solutions is not feasible as it is assumed the Meyerhoff theory which considers the elastic halfspace.


Introduction
The limit load to shallow foundations and the corresponding displacements and material states is one of the most investigated areas in the field of Geomechanics. Kötter and Prandtl [26,45] were the first to investigate the problem providing analytical solutions to sliding line equations and bearing capacity in the context of the infinite halfspace soils. [35,60] investigated the bearing capacity of soils taking into account the friction angle, the soil depth and soil density. Subsequently, a number of papers have been published which investigate the problem [18,29,38,58,66]taking into account the failure mechanism and applying Mohr-Coulomb material yield criterion alongside with linear elastic analysis. In most cases 1D or 2D physical problems have been examined taking into account homogeneous or layered soils: [40,47]. All studies presented were applied to the foundation design regulations by introducing the variables N and S which are for material influence and shape influence, respectively. The material related variables N are three respecting for the 3 addends of the total bearing capacity. N q , N c , N are defined for the influence of possible vertical load in the lateral of the foundation, the cohesion of the soil and the footing dimensions alongside with the total weight of the soil, respectively. In a similar manner the shape variables S q , S c , S are determined.( [31,36,37]) The uncertainty quantification of the bearing capacity of saturated porous media in relation with the input variability such as the Young Modulus or the permeability is investigated with the implementation of the stochastic finite element method [17,44]. In the aforementioned method, the input variability can be simulated by using two alternative methods [8,12,20,23,30,32,33,41,52]. The first method consists of the nodal points to be assumed as random variables and deterministic shape function are implemented for the interpretation of the material spatial distribution [5,14,28,58]. An alternative approach is the implementation of random field series expansion such as the spectral representation or the Karhunen-Loeve expansion or the spatial average method for providing the material input variability [1,4,13,41]. In both approaches, the sampling can be pseudorandom and non-biased or importance sampling methods can be implemented such as the Latin Hypercube Sampling (LHS) [3,53]. Then, for both methods, the standard Monte Carlo simulation can be implemented. In some cases there are studies that implement realizations of non-Gaussian random fields to the input material variables [43].
In the related studies the problems and methodology proposed refer to 1D-2D elastic halfspace theory with material criterion the Mohr-Coulomb yield stress model. In the present work a numerical simulation approach with a modified Cam-Clay material yield model proposed by [24] is adopted, which is an accurate and reliable quantitatively material for clayey soil simulation ( [63]).This material model with the implementation of the finite element method can be used in general 3D loading and geometric circumstances. The excessive computational cost which is required by the Crude Monte Carlo simulation method, can be reduced by efficient computational methods as suggested in [55,56]. For the calculation of failure load, a recurrence relation algorithm which provides accurate results with a relative small amount of trials and with one initial guess trial is implemented. The algorithm is theoretically set, proved and compared with the classic intersection method. The goal of the present paper, which is an extension of previous work by the authors [51], is to quantify the uncertainty of the failure load and displacement alongside with the uncertainty quantification of the failure mechanism in relation to the variability of the input parameters, i.e., the spatial distribution of the material variables and the soil depth. Soil depth is defined as the distance between the upper layer of the soil and the impermeable and rigid bedrock. As the soil depth increases the total stiffness of the soil mass is reduced, leading to smaller values for failure loads. Moreover, with the increase of the soil depth the Meyerhoff spline can be expanded to a more generalized failure mechanism in comparison with localized failure surfaces or punching shear failure( [15,19,61,62]).
The input variability refers to three material variables, namely the compressibility factor , the permeability k and the critical state line inclination c of the soil. The possible spatial distributions of the material variables with respect to the depth Z of the soil domain are the constant, the linear variation and the random field process calculated through the Karhunen-Loeve series expansion with an exponential autocorrelation function. In the constant and linear variation with respect to depth, the truncated normal random variable [7,49] is assumed at the nodal points. It should be noted that LHS importance sampling method is adopted for the construction of the input random vectors. The geometrical dimensions of the 3D soil domain and correlation lengths are parametrically studied and compared with respect to the corresponding solid problem without considering the pore pressure.

Formulation of the dynamic soil-porefluid interaction and the numerical solution
The Biot system of equations governs the physical behavior of porous media static or dynamic loading. If the dynamic load is of low frequency or if static loading problems are analyzed the aforementioned system is reduced leading to a notable decrease of the computational effort. The u-p formulation of Biot system of equations consists of soil-fluid momentum balance with the Darcian flow, the stress-strain law and the boundary conditions. In the present article static forces are applied to the clay soil domain consequently the u-p formulation is adopted. The finite element discretization of the u-p formulation takes the form: where where is the matrix of permeability n units [length] 3 [time]∕[mass] , is the loading vector divided by the total mixture density and Q is a factor influenced by the bulk moduli of fluid and soil skeleton.
are the shape functions of pore pressure in matrix representation.
, , are the standard mass, damping and stiffness matrices of the solid skeleton. Furthermore, , , are the coupling, permeability and saturation matrices,

Definition of plastic and bond strength envelope
The material model,which here is presented in a brief description, is a modified Cam-Clay type model. All the stresses in the following equations of Sect. 3.1 are effective stresses in the solid skeleton. The implemented model is a two surface model, namely the plastic yield envelope (PYE) which defines the elastic region and the bond strength envelope (BSE) which defines the region in which PYE can be located as PYE is always inside BSE. [9,10,22,24,63]. BSE is directly influenced by the structure of the microtiles of the clay. If a stress point lies in BSE, the structure degradation rate of the clay is maximum. Both envelopes have ellipsoidal shape as depicted in Fig. 1. The mathematical representation of BSE is given by the following: where p h is the hydrostatic component of the stress tensor, s is the deviatoric component of the stress tensor, c is the critical state line inclination and a is the half-size of the large diameter of the ellipse Furthermore, PYE is described by the function where p L is the hydrostatic component of the centroid of PYE , is the deviatoric component of the centroid of PYE and is the similarity factor Finally, inside PYE the isotropic poroelasticity theory applies, where the bulk and shear moduli have constant ratio, with the assumption of constant Poisson ratio. The bulk modulus is given by: where is the specific volume of the soil.
The constitutive model is valid for clays, for static and dynamic loads, regardless of the overconsolidation ratio. Clays of friction angle between 17 o -30 o can be simulated with substantial accuracy. This range of friction angles corresponds to the majority of the natural clayey soils. Furthermore, it is proved to be numerically stable because the majority of the equations of the criterion are in a closed form. Furthermore, when highly consolidated clays are investigated it can be easily transformed to take into account a possible cementation of the soil, i.e., to bear tensile stresses and the BSE to move to the left. It should also be noted that the stresses and the strains transformations are made in order to have energy conjugate amounts by applying the numerical transformations used in Von Mises yield criterion. To this context, the following variables are introduced q s = ∶ .
The variable q s stands for the deviatoric stress measure The input material variability can be simulated either assuming that the nodal points are random variables with deterministic shape functions, or by implementing the Karhunen-Loeve series expansion for the construction or a stochastic process [17,23,27,39,42,46,64,65]. In the present study, both approaches are adopted in order to evaluate each approach influence. For the first technique, which is proposed in [64] the random function f is estimated with the usage of shape functions N i by where N 0 is the total number of shape functions and the f i are the values of f in the nodal points which can be random variables following a probability density function (PDF). In the present work, N i are linear functions and the f i follow the truncated normal distribution [2,6,7,49] which has a pdf described by the equation where (X 0 ) and Φ(X 0 ) are the standard normal probability and cumulative distribution function for X 0 , respectively, and A, B, X 0 are the normalized coordinates of the subspace limits and x, respectively. Also, d and the mean value of the normalization refers to the pdf of the random variable before truncation.
For the Karhunen-Loeve implementation, let H 1 ( , ) be a random field of mean ( ) based on a known autocovariance function C h ( , ) = ( ) ( ) ( , ) , where ( , ) is the correlation function and ( ) is the standard deviation of . Therefore, any realization H 1 of the field with M number of eigenfunctions i with corresponding eigenvalues i can be expanded as: where i is a set of random variables of zero mean and covariance function COV ( i , j ) = ij . Finally, for a Gaussian random field, as implemented in the present study, the i functions are a set of standard normal random variables. This type of expansion is the most common because it is strong and robust. In the present work the Karhunen-Loeve expansion is applied by adopting an exponential autocovariance function which has an analytical solution of the Fredholm eigenfunction problem [23,46] where b is the correlation length. If the Fredholm equations cannot be solved analytically, and this occurs when the autocovariance function is more complicated, numerical methods can be applied such as the Galerkin method [20,46]. In the present work, the input material variability is interpreted with both approaches discussed in this section. In the Karhunen-Loeve series expansion, the spatial variability is assumed only in vertical direction.

The Latin hypercube sampling
The selection of the samples of a random variable can be used by crude Monte Carlo simulation, using pseudorandom algorithms to take a large vector of random values, or variance reduction methods can be implemented such as the importance sampling [11] and the Latin Hypercube Sampling (LHS) [3,53]. The implementation of LHS method saves a substantial computational effort in order to estimate mean value and coefficient of variation.
In the Latin Hypercube approach let X be a random vector ( x 1 , x 2 ,...x n ). The n dimensional LHS method is stated as follows: For each random variable x i , the interval [0,1] of the cumulative distribution function (CDF) into N equal subintervals is intersected. Then from each subinterval a random number is chosen and through the inverse CDF a sample of x i is obtained. Once samples for all subintervals and all the random variables are acquired then the x i vectors are randomly permuted and create the vector realization X. With this procedure, it is assured that to each possible row and each possible column in the nXn euclidean space exactly one sample is taken. The visualization of this feature is depicted in Fig. 2.
The advantage of this procedure is that a reduced amount of values in comparison with crude Monte Carlo approach are needed to integrate the PDF of the input and subsequently to estimate the variability of the output. Also, the subintervals in each dimension may not be equal thus taking into account possible asymmetries of the PDF of the input. Finally, this approach can be applied also to cross-correlated variables, i.e., when the correlation matrix is not diagonal thus it has a general use to the uncertainty quantification scientifical field. In the present work the variables are not cross-correlated and the LHS sampling is used to obtain the variables in the three dimensional space of the compressibility factor k, the critical state line inclination c and the permeability k. The above assumption is reasonable because from experimental data available in the literature the cross-correlation of and k can be set to zero due to the coefficient of variation of permeability can be found in some cases much higher than that of volume compressibility, while the correlation between and c can be neglected due to lack of experimental data that provide a correlation relation although indirect empirical expressions for the aforementioned variables have been conducted ( [21,48]).
The material random variables expressed by the PDF g 1 of Sect. 4.1 and the random fields realizations that occur from eq. (9) influence the finite element system of equations (1). The corresponding matrices , , and are changing due to the randomness of the compressibility factor , the critical state line inclination c and the permeability k. The selection of the samples follow the importance sampling method of Latin Hypercube Sampling Method.

An algorithm for the determination of failure load in ramp dynamic load function
In this section, an algorithm for the determination of failure load, when the dynamic loading function is the ramp loading function, is presented. Ramp loading time function is when the load-time relation is linear with zero initial load until time T. The aim of the algorithm is to find the failure load at the time of the end of the ramp loading. Assume the load function of the problem to be as depicted in Fig. 3.
The aim of this algorithm is to define the load factor * which causes failure of the continuum at exactly the time T. An initial guess of 1 is taken leading to an initial time of failure t 1 is obtained. Then, for each new trial n+1 if the load factor n causes failure to the continuum is given by equation Otherwise, if the load factor n is not causing failure, then the maximum no failure factor max−no−failure is obtained from all previous trials implemented for the calculation of n+1 , as follows: It should be noted that in practice this recurrence relation usually converges "by the failure region" which means that only relation (11) is implemented. The difference between n+1 and n is given by: In eq. (13) it is obvious that as n → ∞ , t n → T consequently, the left side of the equation tends to 0, thus the algorithm converges to the desired load factor * .
In comparison with the typical bisection method, where one should guess initial values of failure, safety values 1,fail and 1,no−failure and then calculate the new load factor by the bisection of maximum safety factor and minimum failure factor, this algorithm provides less or equal number of trials for the same initial failure guess and convergence tolerance, as can be seen in Table 1. Moreover, it needs only one initial guess, which in high uncertainty problems is very useful in general for avoiding divergence of the solution. Finally, as proven by the numerical tests presented in Table 1 and in Sect. 6, the difference in failure load and failure displacement between the two algorithms is less than 1%. The absolute percentage difference is computed considering as exact solution the bisection algorithm with initial value of safety stress 100 KPa. The performance in terms of the computational time for the 100 deterministic analyses of a Monte Carlo simulation is also depicted in Table 1. The proposed algorithm can save up to 35% of the time demanded for a Monte Carlo analysis to be performed indicating a notable advantage of the aforementioned recurrence relation.

Description of the problem
The proposed numerical simulation model, which is defined by the set of equations (1) is implemented to porous problems. The geometry of the problem is portrayed in Fig. 5. A uniform vertical load q is applied to the whole soil domain and the ultimate value q fail alongside with the ultimate displacement u fail of point A in Fig. 5 are the monitored output variables. The values for the depth are chosen h = 20 , 40, 50 m. The finite element mesh is with 8 node hexahedral finite elements with linear shape functions for u and p, which provides high numerical accuracy [34,59]. The length in X and Y directions are l x = l y = 4h . The stresses due to geostatic loading are directly imported as initial conditions with the relations v = z , x = y = 0, 85 v as portrayed in stress point L of Fig. 1. The duration of the simulation in all cases is 1 day in order to have quasi static conditions and in each load case a time step of dt=0.001 day is implemented. The other deterministic properties of the soil are given in Table 2.
Here, 0 stands for the initial specific volume of the soil and c stand for the inclination of Isotropic Compression Line for the respective virgin normally consolidated clay and is considered proportional to . The following equations apply to boundary surfaces: x (z = h) = y (z = h) = z (z = h) = , and the lateral boundary surfaces are free of constraints. The input material uncertainty consist of the material variables of the compressibility factor , the critical state line inclination c and the permeability k.
For the compressibility factor , the spatial distribution along the depth is considered as linear denoted as L , or constant denoted as C . For L , depth=0 = 0, 008686 and the ratio R = depth=max depth=0 follows the truncated normal distribution.
The mean value of the ratio is R = 0, 469 and the corresponding CoV= R R is 0,25 , as a consequence z=max,mean = 0, 004074 and the CoV increases with depth. These values are adopted in order for the mean stiffness of the soil to correspond to a shear velocity of 200 m s . In Fig. 4, the linear spatial distribution for is depicted where the realizations of the compressibility factor are depicted for the values of R: R , R + R , R − R . At this point, it is emphasized that that bulk and shear moduli are assumed proportional, since Poisson ratio is assumed constant, consequently is directly associated with the shear velocity. In the case of C , the mean value of is = 0, 004074 and the CoV is 0,25.
For the critical state inclination c the spatial distribution with respect to depth is assumed constant. Two possible cases for the absolute value are considered: the random variable analysis and the deterministic value. In the random variable analysis denoted as c R , the friction angle follows the truncated normal distribution PDF g 1 of Sect. 6. The mean value is = 23 o while the standard deviation is = 2 o and this satisfies the condition that the friction angle to be within an acceptable closed space for physical clays [24].
Subsequently, values of are generated and c is computed from c =  Table 3, associating linear (L) or constant (C) distribution for and deterministic (D) or random variable (R) cases for c. The porous analyses simulated are portrayed in Table 4, combining linear (L), constant (C), and random field (RF) distribution for , deterministic (D), random variable case (R) and random field (RF) distribution for c and random field distribution (RF) for k.
In the random field processes, the mean values are a s s u m e d : mean = 0.008686 , c mean = 0 , 7 3 3 6 a n d k mean = 10 −8 m 3 s Mgr [50,57,67]. The standard deviations are adopted: = 0.25 mean , = 2 o and k = 0.25k mean . The autocorrelation function depicted in 10 is implemented in all stochastic processes. The values of correlation length are b=75 ( k RF� ▽ ) and b=100 ( k RF∞′′ ). The L and C spatial distributions for as well as the random variable distributions for all material variables correspond to a random variable case analysis. In addition, for c a constant deterministic analysis is adopted. The random field (RF) distributions refer to the Karhunen-Loeve series expansion and realizations of the spatial stochastic process are computed through equation H 1 of Sect. 6 with the usage of the autocovariance function 10.
The simulations are static, while the number of Fredholm eigenfunctions taken into account is eight. Failure is defined in the present work as when the first Gaussian Point is having softening behavior (i.e., plastic hardening modulus H0). Each Monte Carlo simulation was applied for 100 samples using the Latin Hypercube Sampling method, which were found sufficient in achieving convergence for the mean value and standard deviation of the monitored displacements as it is depicted in Fig. 6. Here, it should be noted that the cross-correlation in all Research Article SN Applied Sciences (2021) 3:659 | https://doi.org/10.1007/s42452-021-04631-3 material variables is set to zero and as a consequence the correlation matrix is diagonal.

Ultimate failure load and displacement
The results are depicted in Tables 5, 6, 7 and in Figs. 7, 8, 9, 10, 11, 12. In these tables the mean values q fail and u fail of point A in Fig. 5 are presented together with the standard deviation,the coefficient of variation (CoV), the maximum and minimum value of the Monte Carlo simulation.
When the pore pressure is ignored, larger mean failure displacements and smaller CoV are obtained when L is assumed compared to C , while larger mean failure load and corresponding variability are obtained when C is assumed as can be observed from Table 5 for all depths considered. In the numerical tests performed, the largest CoV of the output for q fail is almost half the input variability while the corresponding CoV of the output for u fail is practically the same as the input variability of 0,25 (see -C -c R for h = 20 m). The mean value of u fail in -L -c D for h = 50 m is 48 % larger than the corresponding value for h = 20 m while the mean value for q fail in -C -c D for h = 20 m is 20 % larger than the corresponding value for h = 50 m. Consequently, the critical spatial distribution of for the CoV of the output of both q fail and u fail is the C case. The PDF's of solid analyses are depicted in Figs. 7 and 8 . This behavior in failure displacements can be explained by the fact that in the L case the upper layers of the soil, which are the most compressible, have low CoV of leading to low CoV for both displacements and strains, while in constant distribution along the depth the soil domain tends to be more homogeneous and stiff thus more Gauss Points have larger stiffness leading to larger failure loads.
The CoV of the failure load in porous analyses is slightly influenced by the change in the considered   Table 6. The largest output CoV of q fail , which is found at − C -c R -k RF� ▽ for h = 20 m, is 47 % lower than the input variability, while the corresponding maximum CoV of u fail is 26 % greater than the input variability of 0,25 and is located at − C -c R -k RF∞′′ for h = 20 m. Thus, when considering the pore pressure in the soil domain a reduction of variability of the failure load occurs and in failure displacements when there is constant spatial distribution for a significant variability increase take place as is portrayed in Figs. 9 and 10 . This is in accordance with the results obtained in [17]. This way of behaving is attributed to the fact that when constant spatial variability of a material variable is assumed all Gauss Points are provided with larger uncertainty leading to greater CoV for strains and displacements.
In addition, taking into account that the bulk modulus in porous problems as seen from equation 7, are generally smaller than the respective solid problems, it is concluded that smaller mean values of failure load are expected and smaller variability due to tensile failure of the first Gauss Point as it will be proven numerically in Sect. 6.2.2. and since the tensile stress is set to 0 there are limited acceptable positions for BSE leading to the aforementioned results. The porous analysis with random field representations for all material variables is subsequently performed as a more general case, since it takes into consideration the spatial randomness of the material properties of the soil. In the case of porous random field analyses, the largest CoV of the output for q fail is 3,5 times the input uncertainty, while for u fail the largest output variability is 2,2 times the CoV of the input. The mean values of failure load in porous random field analyses are significantly smaller than the corresponding porous random field analyses with deterministic shape functions for and c, while the mean values for failure displacement are notably greater when the linear distribution over depth is assumed for as is depicted in Table 7 in − RF -c RF -k RF∞′′ of h = 50 m and − L -c D -k RF∞′′ of the same depth. In porous random field analyses          the increase of correlation length decreases the CoV of the output in both q fail and u fail in depth 20 and 40 meters and reverses in depth 50 m.The opposite occurrence is present for the mean failure load. These behavior can be explained by the fact that the types of failure of Gauss points, when the proposed material yield model is applied, are two: the failure from the "wet" side which in Fig. 1 is the stress points laying in the right side of the intersection between BSE and the critical state line and the "dry" side which is the stress points laying in the left side of the same intersection.
Considering the large uncertainty of c and triaxial loading of the Gauss points the failure stress states may be significantly different leading to larger uncertainty of the output. Consequently, the critical spatial distribution for smaller mean failure load, which is the unfavorable situation in the present work, is the Karhunen-Loeve random field representation for all three input variables, which in most cases it also provides the greater output uncertainty.The PDFs of − RF -c RF -k RF analyses are portrayed in Figs. 11 and 12 . The results acquired by the analyses examined depict the effect of the input uncertainty of each material variable in porous failure problems. The compressibility factor , influences the mean value and the variability of both monitored variables. This influence in CoV of the output is more notable when the distribution of with respect to depth is constant. This holds in both porous and nonporous problems. This behavior can be interpreted by the fact that is directly related with the bulk modulus, consequently has a significant influence on the strains, the displacements and the stiffness of the soil domain leading to an influence in the ultimate load.
The permeability k influences to a lesser extent failure load and failure displacement. The spatial variability of k has a slight sensitivity of the CoV of the output in most Monte Carlo simulations for and c with the exception for porous analyses with stochastic processes for all material variables. For the same porous consolidation problem with the same load and deterministic parameters, the output displacement and stress field are not influenced by the permeability since the pore pressures are fully dissipated, leading to the behavior observed by the numerical results.
In conclusion, the variability of critical state line inclination c of the material model appears to have a significant effect on the output variability when the stochastic process calculated from the Karhunen-Loeve series expansion is adopted, while greater mean values are obtained when the c R case is implemented. If a deterministic value for c is chosen the output variability is negligible leading to a conclusion that the influence of this material variable is the most important in comparison with and k. This is explained as c is directly associated with the failure state of the Gauss Point, leading to influence the limit state of the soil mass.
In order to prove that the output displacement and failure load follows the truncated normal or the lognormal distribution, the histograms of three randomly selected Monte Carlo simulations for an output monitored variable with the corresponding distribution fitting are presented in Fig. 13. As proven by the diagrams, the empirical PDF estimated by the histograms can be modeled by the truncated normal PDF g 1 described in Sect. 6 or the lognormal PDF, respectively. The limit values indicated in Tables 5, 6, 7 are accurate values to set the truncation of the PDF. In addition, for proving numerically the Gaussian nature of the output random vectors, the Kolmogorov-Smirnov test is used [16,25,54]. In all output probability density functions, the null hypothesis H 0 , at the 5 % significance level is accepted, as is portrayed in Table 8, where for the Monte Carlo simulation analyses portrayed in Fig. 13 the aforementioned test is implemented. The maximum absolute deviation of the theoretical and empirical cumulative distribution function (CDF) of the Monte Carlo simulations depicted in Fig. 13 are presented and compared to the critical difference for accepting H 0 . In all analyses, the critical value is greater than the maximum absolute difference of the CDFs consequently, the output random vectors follow a Gaussian type distribution. This comes in agreement with the previous research literature [17].

Ultimate stresses-strains and failure mechanism
The results are portrayed in Tables 9, 10, 11. In these tables the mean value, the CoV and the minimum values for the volumetric stresses denoted with p h , the deviatoric stresses denoted with q s , which is the Von Mises stress explained in Sect. 3.1. In addition, the mean value with the CoV of the volumetric strains denoted with vol , deviatoric strains denoted with e are also depicted in Tables9, 10, 11. Moreover, the probability of first Gauss Point failure for all simulations is presented. Here, it should be noted that due to symmetry of the displacement field only the point with smaller coordinates is written. It is also pointed out that the volumetric strain at failure is tensile. In solid analyses, larger CoV of the output and smaller minimum value of failure stress is obtained when c R case is adopted as it can be seen from Table 9. In most cases, the aforementioned distribution also provides greater mean values. The same spatial distribution assumption for critical state line inclination provides larger variability of the strains at failure, which on average are on the proximity of 3-4 ‰. In general, from the results of the Monte Carlo simulations it can be concluded that the percentage plastic deviatoric strains are larger leading to the conclusion that the distortional failure is critical. Finally, in most cases the Gauss Point (2,11 , 2,11 , 12,11) is the most probable failure point with the exception of the c distribution at h = 40 m.
In porous analyses with deterministic shape functions for and c larger stress variability is obtained at h = 20 m with deterministic value for c and at greater depths with c R case. The change of the correlation length influences notably the output variability at h = 20 m as is portrayed in Table 10. The maximum output variability in failure strains occurs at h = 50 m and correlation length for k b = 75 m. From the Monte Carlo simulations implemented, it can be seen that when the c -c R distributions are assumed the volumetric failure is critical, while when c -c R distributions are adopted the deviatoric failure occurs. In conclusion, the most probable failure point is the (2,11 , 2,11 , 12,11) for h = 20 m, the (2,11 , 2,11 , 32,11) for h = 40 m and the (2,11 , 2,11 , 2,11) for h = 50 m and linear distribution for , while (2,11 , 2,11 , 42,11) is the unfavorable Gauss Point    for h = 50 m and constant distribution along the depth for the compressibility factor. In porous analyses with random field representation for , c, k larger output CoV in stresses is obtained in general when b = 75 m except from the deviatoric stresses in h = 40 and 50 m as depicted in Table 11. The volumetric strains and the deviatoric strains have in all depths critical correlation length b critical,vol =75 m and b critical,dev =100 m, respectively. On average, considering the numerical results obtained by the aforementioned analyses, at h = 20 and 50 m the volumetric failure is critical, while at depth 40 m the distortional failure occurs. Finally, the most probable failure points in 20 m are the (2,11 , 2,11 , 2,11) for b = 75 m and (2,11 , 2,11 , 12,11) for b = 100 m, while in depth 50 m is the (2,11 , 2,11 , 32,11) for both correlation lengths. If the soil depth is 40 meters, there are many points with equally high probability of failure thus leading to a more uncertain position of the failure mechanism.

Conclusions
In this paper, the uncertainty quantification of the failure of clayey soils taking into account the pore pressure-soil interaction with the implementation of the stochastic finite element method is presented. The aim of this work is to introduce accurate quantitative results on the failure load of 3D clayey soil domain in relation with the input variability of soil material variables. The proposed numerical simulation model holds for every possible assumptions for the geometry, the loading properties and the material  Table 9 Monte Carlo results for the stresses (KPa), the strains and the probability of the first Gauss point failure for non-porous medium                      distribution of the soil domain. In this context, a detailed finite element simulation, alongside with a complicated material constitutive model is used. It is evident that similar computational study can be conducted for shallow foundations and it is a future research work that will be performed by the authors. The numerical results obtained indicate that the output failure load and displacement follow a Gaussian random distribution despite the material nonlinearity, which in failure phenomena is excessive. The randomness of material poroelasticity plays an important role in the output CoV for failure load and failure displacement, especially when it has a Karhunen-Loeve random field representation along the depth of the soil domain. The same thing applies for the critical state line inclination variable c. When the random Karhunen-Loeve field is assumed for all stochastic material variables in discussion the CoV of the output exceeds the corresponding variability of the input in half of the Monte Carlo simulations. This variability increase varies from 30% to 3,5 times larger. The porous variable of permeability influences to a lesser extent the uncertainty of the output material variables. In porous media problems the failure load is smaller compared to the corresponding solid problems.
The random field representations for , c and k provide the maximum variability for failure stresses and strains. In porous random field analyses in depths h = 20 , 50 m the volumetric failure occurs, while when the soil domain is 40 m deep the distortional failure is critical. When deterministic shape functions for and c are implemented, the c -c R distributions provide with larger probability volumetric failure of Gauss Points, while L -c D distributions are adopted in order to obtain with larger probability deviatoric failure. Finally, in most cases (2,11 , 2,11 , 12,11) is the critical Gauss Point of failure thus it can be the starting point of the failure mechanism through the Meyerhoff spline.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creat iveco mmons. org/ licen ses/ by/4. 0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
Funding As stated in the Acknowledgments section.
Availability of data and material Upon communication with the corresponding author.

Conflicts of interest As stated in the Acknowledgments section.
Code availability Open source code MSolve in programming Language C# . Information to the link http:// mgroup. ntua. gr/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.