Uncertainty Quantiﬁcation of Failure of Shallow Foundation on Clayey Soils with a Modiﬁed Cam-Clay Yield Criterion and Stochastic FEM

: In this article, a quantitative numerical study of the random distribution of the soil material parameters to the probability density functions of the failure load and failure displacements of a shallow foundation is presented. A modiﬁed Cam-Clay yield function is used for this scope into a stochastic ﬁnite element numerical formulation. Several hypotheses for the random distribution of the compressibility factor κ , of the material constitutive relation, the critical state line inclination c of the soil, as well as of the permeability k of the continuum, have been tested and assessed with Monte Carlo simulation accelerated with Latin hypercube sampling. It is validated that both failure load and failure displacements follow Gaussian normal distribution despite the non-linear behaviour of the soil. Furthermore, as the soil depth increases, the mean value of failure load decreases and the failure displacement increases. The failure mechanism of clays can be determined with accuracy using this numerical implementation, without the restrictions imposed by analytical solutions, taking into consideration the eccentricity of the load in combination with non-linear constitutive relations.


Introduction
The ultimate load of shallow footing settlements and the corresponding material states and displacements is one of the most researched areas in the field of Geomechanics. From the seminal work of Terzaghi [1], who investigated the bearing capacity of soils, taking into account the friction angle, soil depth and the density, a number of papers followed investigating this problem. Among recent ones are [2][3][4][5][6], taking into account the failure mechanism and applying Mohr-Coulomb material yield criterion alongside linear elastic analysis. In most cases, 1D or 2D continua have been examined along with homogeneous or layered soils: [7,8]. All studies presented were applied to the foundation design regulations by introducing the friction variables N and shape-influence variables S. The three friction variables N represent the three factors contributing to the total bearing capacity. N q , N c and N γ are defined for the contribution of possible vertical load in the lateral of the foundation, the cohesion of the soil and the settlement dimensions, along with the total weight of the soil, respectively. Correspondingly, the shape variables S q , S c , S γ are interpreted [9][10][11].
The uncertainty quantification of the ultimate load of shallow foundations of saturated porous media, with respect to the input uncertainty, such as the Young Modulus or the permeability, has been researched with the adoption of the stochastic finite element method, where the input variable and its uncertainty can be interpreted by using two alternative methods. The first method supposes the properties of nodal points as random variables and the shape functions implemented for the interpolation of the material spatial distribution as deterministic relations [4,[12][13][14][15][16]. A different approach is the adoption of random field representation, such as the spectral representation or the Karhunen-Loeve expansion, or the spatial average method, to provide the material input variability [17][18][19][20][21][22][23][24]. In both simulation methods, the sampling can be either non-biased through pseudorandom algorithms, or an importance sampling method could be implemented, such as the Latin Hypercube Sampling (LHS) [25,26]. Then, the standard Monte Carlo simulation can be performed.
Monte Carlo simulations can be applied to other scientific areas of geotechnics. The generality of this method has led to implementation in many different practical applications. Some of them are briefly presented alongside the corresponding scientific publications. For seismic excitations in slopes, as well as for their stability, the computation of the safety factor in relation to the slope angle, the cohesion and the friction angle can be estimated through Monte Carlo analyses [27]. In the ground, vibrations occurring from blasting the risk estimation can be assessed through stochastic procedures and Monte Carlo algorithms [28]. Moreover, in mine engineering and hard rock pillars, a stability analysis can be conducted through stochastic finite element method and reliability assessment can be obtained [29]. Finally, a prediction of flyrock distance due to blasting can be obtained, assuming a random distribution of the material of the rock domain and implementing a Monte Carlo simulation technique [30].
Previous studies refer to 1D-2D elastic halfspace theory with material criterion as the Mohr-Coulomb yield stress model. In the present paper, a numerical simulation perspective with a Modified Cam-Clay material yield model proposed by [31] is implemented, which is a precise and quantitatively valid material for clayey soil simulation [32]. This material yield function with the adoption of a finite element model can represent the real load and displacement field in every possible 3D-loading condition. The immoderate computational cost, which is needed by the crude Monte Carlo simulation method, can be alleviated with efficient computational schemes, as suggested in [33,34]. For the calculation of failure load, an improved version of a recurrence-relation algorithm proposed by the authors [35], which yields valid results with a small number of trials and with one initial guess trial, is presented. The algorithm is theoretically set, proved and compared with the classic bisection method. The goal of the present paper, which is an extension of previous work by [35], is to quantify numerically, in the case of shallow foundations, the variability of the failure load and failure displacement alongside the uncertainty quantification of the failure spline in reference to the uncertainty of input parameters such the spatial distribution of the material variables and the soil depth. The aforementioned results provide a reliability analysis for the variability of the failure load, failure displacements and the corresponding mechanism of the Meyerhoff spline. In addition, this study offers a numerical investigation algorithm for fastening the procedure of failure-load estimation.
The numerical investigation consists of stochastic analyses of footing settlement failure in clayey soils. A rectangular soil domain subjected to the equivalent forces of a shallow foundation under single eccentricity is assumed, and three material variables are considered uncertain: videlicet, the compressibility factor κ; the permeability factor k; and the critical state line inclination c of the soil. The spatial distributions of the material parameters in relation to the depth z of the soil domain are the linear, the constant variation and the random field process, constructed from the Karhunen-Loeve series expansion, considering an exponential auto-correlation function. In the constant and linear variation with respect to depth, the truncated normal random variable [36,37] is presumed at the nodal points. The LHS importance sampling method is adopted for obtaining the input random vectors. The eccentricity of the footing settlement is parametrised with four possible values and correlation lengths are parametrised with three possible values and contrasted to the analogous solid problem in which the pore pressure is neglected.

Formulation of the Dynamic Soil-Pore-Fluid Interaction and the Numerical Solution
The response of porous media under static and dynamic loading conditions can be predicted with the Biot mathematical system of equations. When a low-frequency load is applied or static loading is considered, the formulation is simplified, leading to a considerable alleviation of the complication of the problem and, subsequently, of the computational time needed. The u-p formulation of Biot system of equations comprises the soil-fluid momentum balance alongside the Darcian flow, the stress-strain law and the boundary conditions. In the present work, the u-p formulation is adopted, since static forces are implemented to the clay soil domain.
The finite element discretisation of the u-p formulation takes the form: where k is the permeability matrix, b 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. N p are the shape functions of pore pressure in matrix representation. M st , C st , K st are the standard mass, damping and stiffness matrices of the solid skeleton. Furthermore, Q c , H, S are the coupling, permeability and saturation matrices, respectively. Finally, f eq corresponds to the equivalent forces due to the external loading. Numerical schemes, such as the Newmark direct integration method, are implemented to obtain the solution to the problem.

Plastic Yield Envelope and Structure Strength Envelope Mathematical Representation
The material yield function used for describing the material stresses in the solid skeleton in this article is a modified Cam-clay-type model. This is a two-surface model: the plastic yield envelope (PYE) for the definition of the elastic area and the structure-strength envelope or Bond-Strength Envelope (BSE), which defines the limit states that PYE can have [31,32,[38][39][40]. BSE is mainly affected by the structure of the cohesive soil. If a stress point lies in the BSE boundary, the structure degradation rate of the clayey soil is the greatest. Both envelopes are assumed as ellipses and are portrayed in Figure 1.
The mathematical representation of an envelope is the following: In Equation (5), the general stress point σ has a hydrostatic component p h and a deviatoric component s, while the centre of the ellipse L has a hydrostatic component p L and a deviatoric component s L . Moreover, a is the halfsize of the large diameter of BSE, and a similarity factor, ξ, is introduced. When the deviatoric part of the centre of the ellipse is not zero and the hydrostatic part is not equal to a, then the generalised envelope coincides with the plastic yield envelope of the stress point.
However, if s L = 0 , p L = a then ξ = 1 and, consequently, the Bond-Strength Envelope is obtained: For the simulation of the elastic behaviour the integration point is assumed poroelastic. The bulk modulus, proportional to the shear modulus as a consequence of a constant Poisson ratio, is as follows: ν denotes the specific volume of the soil. This material law for the majority of natural cohesive soils is substantially accurate, especially when the friction angle is in the range of 17 • and 30 • . The reliability of the simulation of the constitutive behaviour holds for all possible loadings. In addition, since the equations of the material-constitutive model are in closed form, the derived numerical scheme is stable. Moreover, when large values of overconsolidation ratio (OCR) of a cohesive soil are simulated, this yield function can be easily reconstructed to take possible tensile stresses into consideration. The stresses and strains' algebraic transformations are implemented in order to have energy conjugate amounts by adopting the numerical transformations used in von Mises yield criterion. In this context, the following variables are introduced. q = √ 3/2s : s , e = 2/3 dev : dev (9) The variable e is used to denote deviatoric strain measure, whilst q is used to denote the deviatoric stress measure (Von Mises stress).

The Karhunen-Loeve Series Expansion and the Truncated Normal Distribution
The material uncertainty of the input can be interpreted either by considering the nodal points as random variables with deterministic shape functions, or by performing the Karhunen-Loeve series expansion for the computation of a realisation of a stochastic process [41][42][43][44][45][46][47]. In this work, both approaches are adopted in order to assess the perfor-mance of each simulation method. For the first approach, which is proposed in [41], the random function f 0 is estimated with the use of shape functions N i : (10) where N 0 is the total number of shape functions and the f 0,i are the nodal point values, 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 [36,37,48,49], which has a PDF described by the equation In Equation (11), φ(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 normalised coordinates of the subspace limits and x, respectively. Additionally, σ d and the mean value of the normalisation refers to the PDF of the random variable before truncation.
For the Karhunen-Loeve series expansion, let H 1 (x, ω) be a random field of mean µ(x) based on a known auto-covariance function where ρ(x 1 , x 2 ) is the correlation function and σ(x 1 ) is the standard deviation of x 1 . Any realisation 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 E(ξ 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 widely used due to its high efficiency. In the present paper the Karhunen-Loeve expansion is applied by adopting an exponential auto-covariance function which has an analytical solution of the Fredholm eigenfunction problem where b is the correlation length. If the Fredholm equations cannot be solved analytically, and this occurs when the auto-covariance function is more complicated, numerical methods can be applied [18,44].

The Latin Hypercube Sampling
The samples of the random variables can be chosen by implementing pseudorandom unbiased algorithms to generate a large random vector, or by variance-reduction methods such as the importance sampling [50] and the Latin Hypercube Sampling (LHS) [25,26]. The implementation of the LHS method preserves a notable amount of computational time in order to estimate the statistical moments of the output random variables.
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) is intersected into N equal subintervals. 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 realisation X.
With this procedure, it is assured that at each possible row and each possible column in the (n × n) Euclidean space, exactly one sample is taken.
The ascendancy of this scheme is that a smaller amount of values in contrast with the crude Monte Carlo simulation are needed to integrate the PDF of the input and, consequently, to estimate the uncertainty of the output. In addition, the subintervals in each dimension may not be equal, thus taking into account possible asymmetries of the PDF of the input. In conclusion, this methodology can also be implemented to cross correlated variables, i.e., when the correlation matrix is not diagonal; thus, it has a general use to the uncertainty quantification literature. In the present study, 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 material random variables expressed by the PDF g tr of Section 4.1 and the random field realisations experienced in Equation (13) influence the finite element system of Equation (1). The corresponding matrices C, K and F are changing due to the randomness of the compressibility factor κ, as well as the critical state line inclination c and the permeability k. The selection of the samples follows the importance-sampling method of the LHS technique.

An Improvement of a Proposed Algorithm for the Determination of Failure Load in Ramp Dynamic Load Function
In this section, an improvement of an algorithm proposed by the authors [35] for the determination of failure load, when the dynamic loading function is the ramp loading function, is presented. The aim of the algorithm is to find the failure load at exactly the end of the ramp loading, as indicated in Figure 2. 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, and 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, it is given by the equation where a and b are weight coefficients chosen in order to improve the convergence speed. If the load factor λ n is not causing failure, then the maximum no-failure factor λ max−no− f ailure 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 (15) is implemented. The difference between λ n+1 and λ n is given by: In Equation (17), it is obvious that as n → ∞, t n → T. Consequently, the left side of the equation tends to zero, and subsequently, the algorithm converges to the desired load factor λ * if and only if In Table 1, an investigation of the weight coefficient's influence on the convergence speed is presented, with reference to the typical bisection method. In comparison with the typical bisection method, where one should guess initial values of failure, safety values λ 1, f ail and λ 1,no− f ailure , and then calculate the new load factor by the bisection of maximum safety factor and minimum failure factor, this algorithm provides a lower 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 Section 6, the difference in failure load and failure displacement between the two algorithms is less than 1%. The absolute percentage difference is computed considering the bisection algorithm with initial value of safety stress 1000 kPa as an exact solution. 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 from almost 5% to almost 55% of the time demanded for a Monte Carlo analysis to be performed, in the case of a = 1 and b = 0, indicating a notable advantage of the aforementioned recurrence relation in all the family of relations, indicated by Equation (15). Table 1. Comparison of the proposed algorithm with the bisection method for the calculation of the failure load λ * and failure displacement.

Description of the Problem
The proposed numerical simulation model is implemented in porous problems as portrayed in Figure 3 and is defined by the set of Equation (1). The monitored output variables are the maximum settlement stress, considering linear variation over the surface of the foundation, the normal force of the footing settlement, the maximum displacement of the settlement-which is considered as the mean output value of the points A and B of Figure 3-and the minimum displacement of the settlement, which is considered as the mean output value of the points C and D of Figure 3. The loading conditions consist of the nodal values q 1 at points C and D and q 2 at points A and B and represent, at each loading case, the equivalent forces of a linear distributed loading case for a specific eccentricity of the footing settlement. The eccentricity e = M N is considered with four possible values, namely 0, h 12 , h 6 , h 3 . The finite element mesh consists of eight-node hexahedral finite elements with linear shape functions for u and p, which results in an acceptable numerical accuracy [51,52]. The lengths of the soil domain in X, Y and Z directions are l x = 5 m, l y = 5 m and l z = 4 m, respectively. The stresses due to geostatic loading are directly imported as initial conditions with the relations σ v = γz, σ x = σ y = 100 kPa associated with stress point L of Figure 1. The duration of the simulation in all cases is one day in order to have quasi static conditions, and at each load case, a time step of dt = 0.001 d is implemented. The other deterministic properties of the soil are given in Table 2. In Table 2, ν 0 stands for the initial specific volume of the soil and λ for the inclination of isotropic compression line for the respective virgin, normally consolidated clay, which is considered proportional to κ. The following equations apply to boundary surfaces: u x (z = h) = u y (z = h) = u z (z = h) = 0, and the lateral boundary surfaces are free of constraints. The input material uncertainty consists of the material variables, the compressibility factor κ, the critical state line inclination c and the permeability factor k.
For the compressibility factor κ, the spatial distribution along the depth is considered either linear (κ L ) or constant (κ C ). In the κ L case, κ z=0 = 0.008686 and the ratio R follows the truncated normal distribution with R = κ z=max κ z=0 . The mean value of the ratio is µ R = 0.469 and the corresponding CoV is 0.25; as a consequence, κ z=max,mean = 0.004074.
These values are selected in order for the mean stiffness of the soil to correspond to a shear velocity of 200 m s . In Figure 4 the spatial distribution of κ L is depicted. It should be mentioned that the bulk and the shear moduli are assumed proportional, since the Poisson ratio is assumed to be constant and, consequently, κ is directly associated with the shear velocity. When κ is considered constant, the mean value of κ is κ µ = 0.004074 and the CoV is 0.25. For the critical state inclination c, the constant spatial distribution with respect to depth is adopted. Two possible hypotheses for the absolute value are assumed. If a random variable case is adopted c R , the friction angle φ 0 follows the truncated normal distribution PDF with the mean value of µ φ = 23 • and standard deviation of σ φ = 2 • , which provides values for φ which are in an acceptable closed space for natural clays [31]. Values of φ 0 are generated through the generation of samples from the standard normal distribution with the Latin Hypercube Sampling and then are transformed into the truncated normal PDF g tr . Subsequently, c is computed from c = 2 3 6sin(φ 0 ) 3−sin(φ) . Alternatively, for the deterministic case c D , c = 0.7336 for friction angle µ φ = 23 • .
For the permeability k, the constant spatial distribution with respect to depth is adopted. The absolute value may be considered in two approaches. If a random variable case is adopted k R , the mean value is µ k = 10 −8 and the CoV is CoV k = 0.25. On the other hand, for the deterministic case k D , k = 10 −8 .
Two types of analyses are examined: the solid analyses, where the pore pressure is neglected, and the porous analyses, where the water flow is taken into account. The solid analyses performed, specified with (S), are depicted in Table 3, incorporating linear (L) or constant (C) distribution for κ and deterministic (D) or random variable (R) cases for c. The corresponding abbreviated porous analyses performed are portrayed in Table 4, incorporating linear (L), constant (C) and random field (RF) distribution for κ; deterministic (D), random variable case (R) and random field (RF) distribution for c; and deterministic (D), random variable case (R) and random field distribution (RF) for k.  Horizontal axis stands for the compressibility factor κ and the vertical axis for the depth of the soil point in meters. The spatial distributions for the mean value µ of the ratio R = κ z=max κ z=0 and the corresponding distributions for µ + σ and µ − σ. Table 3. Non-porous (solid) analyses performed.
In the random field processes the mean values are considered: κ mean = 0.008686, c mean = 0.7336 and k mean = 10 −8 m 3 s Mgr following [53][54][55]. The standard deviations implemented are: σ κ = 0.25κ mean , σ φ = 2 • and σ k = 0.25k mean . The auto-correlation function of Equation (14) is adopted in all stochastic processes. The values of correlation length are b = 2 m (k RF−2 ), b = 4 m (k RF−4 ) and b = 8 m (k RF−8 ). The spatial distributions for κ, κ L and κ C , as well as the random variable distributions for all material variables correspond to a random variable case analysis. For c, a constant deterministic analysis is adopted. The random field (RF) distributions refer to the Karhunen-Loeve series expansion and realisations of the spatial stochastic process, and are computed through equation H 1 of Section 4.1 with the use of the auto-covariance function given by Equation (14).
The simulations are static, whilst the number of Fredholm eigenfunctions taken into consideration is eight. Failure is defined when the first Gaussian Point exhibits softening behaviour (i.e., plastic-hardening modulus H < 0). Each Monte Carlo simulation was applied for 100 samples, using the Latin Hypercube Sampling method, which was found sufficient in achieving convergence for the mean value and standard deviation of the monitored displacements, as is portrayed in Figures 5 and 6. It should be noted that the cross correlation in all material variables is neglected, and consequently, the correlation matrix has a diagonal form.

Ultimate Failure Load and Displacement
The results for the statistics of failure load and failure displacement are depicted in Tables A1-A3 and in Figures A1-A6. In these tables, the mean values (µ ν ) of the maximum settlement stress, the normal force of the shallow foundation, the maximum and the minimum displacements and the rotation of the footing are depicted alongside the coefficient of variation (CoV) (σ δ ), as well as the maximum (M) and minimum µ values obtained from the Monte Carlo simulation. It should be noted that in Figures A1-A6, the PDFs of the normal force of the shallow foundation and the maximum settlement displacement are portrayed, while for the maximum footing stress, the minimum displacement and the rotation of the foundation of the corresponding figures are in the supplementary materials section, which is distinct from the main manuscript file.
As can be distinguished from Table A1, when the pore pressure is neglected, greater mean failure displacements and lower CoV are obtained when κ L is assumed compared to κ C , whilst greater mean failure stress and force and corresponding uncertainty are obtained when κ C is assumed for all eccentricities considered, except for eccentricity h 6 and random variable case for c. In the numerical simulations performed, the largest CoV of the output for failure stresses and forces is 40% of the input uncertainty, while the maximum CoV of the output for the displacements and the rotation is almost the same as the input variability of 0.25 (see the variable u min in the S-κ C -c R analysis for e = h 3 ). The mean value of u max in S-κ L -c D for e = 0 is 13% larger than the corresponding value for e = h 3 , while the mean value for φ, in S-κ L -c D for e = h 3 is 3.3 times greater than the corresponding value for e = h 12 . Consequently, the critical spatial distribution of κ for the CoV of the output of monitored stresses, forces, displacements and rotation is the κ C assumption. The PDFs of solid analyses are portrayed in Figures A1 and A2, where it can be concluded that for failure stresses, the mean value is increased with the increase in the eccentricity, whilst in normal forces, the opposite behaviour occurs. For failure displacements, the increase in the eccentricity causes a decrease in the mean value, and this also applies to the settlement rotation. This behaviour in the displacements and the rotation of the foundation can be attributed to the fact that in the κ L assumption, the upper layers of the soil, which are the most compressible, have smaller CoV of κ, causing lower CoV for displacements and strains. The assumption of constant distribution along the depth of the soil mass tends to be more homogeneous and stiff, and subsequently, the Gauss Points have greater stiffness, and greater failure stress and forces occur.
The CoV of the failure stresses and forces in porous analyses is less influenced by the change in the eccentricity, unlike the output variability of the failure displacements and rotation, as it is portrayed in Table A2. The largest output CoV of maximum footing stress, which is found at P-κ C -c R -k D for e = h 3 , is 44% lower than the input uncertainty, whilst the corresponding maximum CoV of u max is 2.6 times larger than the input variability of 0.25 and is located at P-κ C -c D -k R for e = h 12 . Therefore, when taking into account the pore pressure in the soil domain, a reduction in the variability of the failure stresses and forces occurs, whilst in failure displacements, when there is constant spatial distribution for κ, a significant variability increase takes place, as is indicated in Figures A3 and A4. For the rotation of the footing settlement, when the constant distribution of the compressibility factor is assumed, there is a slight divergence from the uncertainty of the input. This behaviour is explained from the fact that when constant spatial variability of a material variable is assumed, all Gauss Points are associated with larger uncertainty, leading to greater CoV for strains and displacements. Considering that the bulk modulus of Equation (8) in porous problems is generally smaller than the corresponding solid problems, it is deduced that smaller mean values of failure load and smaller variability are expected, due to tensile failure of the first Gauss point. This will be demonstrated numerically in Section 6.2.2.
Porous analysis with all the material variables following a stochastic process is performed as a more general case, since it takes into account the spatial uncertainty of the material properties of the soil. In the case of porous random field analyses, the largest CoV of the output for maximum stresses is 36% lower in relation to the input variability, and similar conclusions apply for the maximum foundation displacement u max . For the rotation of the settlement φ, the largest uncertainty of the input is very close to the variability of the input of 0.25, and is found when the eccentricity is h 3 and b = 4 m. The mean values of failure maximum stresses in porous random field analyses are slightly smaller (divergence up to 17%) than the corresponding porous random field analyses with deterministic shape functions for κ, c, k whilst the mean values for maximum failure displacements are greater than the respective analyses with linear distribution for κ, with divergence in the vicinity of 20%, as portrayed in Table A3. In porous random field analyses the increase in correlation length decreases the CoV of the output in failure stresses, forces and displacements in all eccentricities, while for the rotation of the foundation, the maximum uncertainty occurs when b = 4 m. The unfavourable situation arises when the normal force is small thus the critical spatial distribution is considered the one with the smallest mean value as output. Consequently, the critical spatial distribution for the mean value of N is the random field representation for all three material input variables. The PDFs of P-κ RF -c RF -k RF analyses are portrayed in Figures A5 and A6.
The results obtained provide qualitative indications as well as quantitative measurements of 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 all monitored variables. This influence on the 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, which can be attributed to the fact that κ is directly related to the bulk modulus, and consequently has a significant influence on the strains, the displacements and the stiffness of the soil domain, leading to a greater influence of the ultimate load.
The permeability k affects failure stresses, forces and failure displacements and rotations to a lesser extent. The spatial variability of k has a small influence on the CoV of the output in most Monte Carlo simulations for κ and c. 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, as observed by the numerical results.
In conclusion, the variability of critical state line inclination c of the material model appears to have the most notable effect on the output variability for maximum footing settlement and forces when the constant distribution and random variable case is adopted, and greater mean values are obtained when the c R case is implemented. Similar conclusions can be made for the displacements and the rotations of the footing settlement. When a deterministic value for c is chosen, the output variability of the maximum stresses and forces is negligible, leading to the conclusion that the influence of this material variable is the most important in comparison with κ and k. This is explained by the fact that c is directly associated with the failure state of the Gauss point, which influences the limit state of the soil mass.
In order to demonstrate that the monitored output variables follow the truncated normal or the log-normal distribution, the histograms of three randomly selected Monte Carlo simulations for an output monitored variable with the corresponding distribution fitting are presented in Figure A7. As shown in the diagrams, the empirical PDF estimated by the histograms can be modelled by the truncated normal PDF g tr described in Section 6 or the log-normal PDF, respectively. The limit values M and µ indicated in Tables A1-A3 are valid values to set the truncation of the PDF.
In order to numerically prove the nature of the output random vectors, the Kolmogorov -Smirnov test is used [56][57][58]. In all output probability density functions, the null hypothesis H 0 at the 5% significance level is accepted, as shown in Table A4, whereas for the Monte Carlo simulation analyses depicted in Figure A7, the aforementioned test is implemented. The maximum absolute deviation of the theoretical and empirical cumulative distribution function (CDF) computed by the Monte Carlo simulations of Figure A7 are compared to the critical difference for accepting H 0 . It can be seen that 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 [43].

Ultimate Stresses-Strains and Failure Mechanism
The results for the ultimate stress-strain relations and failure mechanisms are presented in tables located in the supplementary materials section. In these tables, the following parameters are considered: mean values, CoV and minimum values for the volumetric stresses, denoted as p vol ; and deviatoric von Mises stresses, denoted as q dev , explained in Section 3.1. Furthermore, the mean value with the CoV of the volumetric strains, denoted as e vol ; and the deviatoric strains, denoted as e dev , are also included in the aforementioned tables, while the probability of first Gauss Point failure for all simulations is also presented. At this point, it is pointed out that the volumetric strain at failure is tensile.
As can be seen in solid analyses, larger variability of the output and smaller minimum value of failure stress is obtained when c R case is considered. In most cases, this distribution also provides larger mean values. The same spatial distribution for critical state line inclination provides larger uncertainty of the strains at failure in almost all eccentricities, which on average are in the proximity of 12‰ for the linear distribution assumption for κ and in the proximity of 5‰ for the constant distribution for the compressibility factor. In general, from the results of the Monte Carlo simulations, it can be seen that the percentage plastic deviatoric strains are greater; consequently the distortional failure is critical. Finally, in all cases, the Gauss point (3.21, 2.21, 3.79) is the most probable failure point when the eccentricity is not zero.
In porous analyses with deterministic shape functions for κ, c, k larger stress variability is obtained at e=0 with deterministic value for c, and when the eccentricity is not zero, the c R case is critical. The change in the permeability notably influences the output uncertainty when the κ C -c R is assumed for eccentricities h 12 , h 6 . The largest output uncertainty in failure strains occurs at e = 0, κ C -c D -k D analysis. From the results acquired with Monte Carlo analyses, it can be concluded that when the eccentricities are small (e = h 12 ), the deviatoric failure is critical, while when larger eccentricities occur, in general the linear distribution for κ provides deviatoric failure and the constant distribution for the compressibility factor provides volumetric failure. In conclusion, the most probable failure Gauss point is the (3.21, 2.79, 3.79) for e = 0 and h 12 . The point (3.21, 2.79, 3.79) is also critical for e = h 6 except from the linear distribution assumption for κ and random variable case for c, whereas the critical Gauss point is the (3.21, 2.21, 3.79). Finally, (3.21, 2.21, 3.79) is the most probable failure point for e = h 3 for all possible assumptions that do not implement random field processes for the material variables in discussion.
In porous analyses with random field representation for κ, c, k, greater output CoV for stresses is obtained in general for b = 2 m. The volumetric strains and the deviatoric strains have, in all eccentricities, critical correlation length b = 4 m. Taking into account the numerical results acquired by the analyses performed, when the eccentricity is greater than zero, the distortional failure occurs. Finally, the most probable failure Gauss point in small eccentricities is the (3.21, 2.79, 3.79) whilst for larger eccentricities, (3.21, 2.79, 3.79) is the critical failure mechanism point for all correlation lengths assumed. As the eccentricity equals to the one-third of the height of the footing settlement, the probability of the most probable point tends to 100%.

Conclusions
In this work, the uncertainty quantification of the failure of shallow foundations on clayey soils is presented, taking into consideration the pore pressure-soil interaction with the implementation of the stochastic finite element method. The aim of this work is to provide a numerical tool for obtaining accurate quantitative results regarding the failure stresses, forces, displacements and rotations of the footing settlement on clayey soil domain in relation to the input uncertainty of soil material variables. The proposed numerical simulation model is valid for every possible assumption for the geometry, the loading properties and the material distribution of the soil domain. In this context, a detailed finite element simulation alongside an accurate and reliable material constitutive model is implemented. The practical application of this article consists of the numerical implementation of the stochastic finite element method through Monte Carlo analysis to a spatially random cohesive soil domain that is subjected to a shallow foundation load that causes failure. The cohesive soil domain has material constitutive model of a Modified Cam-Clay yield function and the output random vectors of failure load, failure displacement and the integration Gauss point that the failure mechanism initially occurs are obtained and assessed.
The numerical results acquired portray that the monitored output variables of maximum stresses, forces, displacements and rotation follow a Gaussian random distribution despite the excessive material non-linearity, which occurs in failure situations. The randomness of material poroelasticity plays an important role in the output uncertainty for failure stresses, forces and failure displacements, especially when it has a constant distribution along the depth of the soil domain. The same observation holds for the critical state line inclination variable c. When the constant distribution for κ is assumed, the CoV of the output maximum displacement exceeds the corresponding variability of the input with an increase which is up to 2.6 times larger. The porous variable of permeability influences the uncertainty of the output material variables to a lesser extent. 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 mean failure strains, especially when the correlation length is smaller. In porous random field analyses for all eccentricities, the distortional failure is critical. When deterministic shape functions for the material variables are considered, for smaller eccentricities the deviatoric failure is the critical. For larger eccentricities of the footing settlement, if the compressibility factor is assumed linear over the depth of the soil domain, the volumetric failure is critical, whilst for a constant distribution for κ, the deviatoric failure occurs. Finally, in most cases (3.21, 2.79, 3.79) is the critical Gauss point failure onset, except from the eccentricity h 3 , whereas the most probable Gauss point for failure is the (3.21, 2.21, 3.79), which can be considered as the starting point of the failure mechanism through the Meyerhoff spline.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/geotechnics2020016/s1, Figure S1: PDFs of the maximum settlement stresses in kPa for non porous medium. κL, κC stands for the linear and constant distribution for κ respectively and c D , c R stands for the deterministic and random distribution for c respectively. Subplot (a), (b), (c), (d) stand for eccentricity 0, h 12 , h 6 , h 3 , respectively; Figure S2. PDFs of the minimum displacement in m for non porous medium. κL, κC stands for the linear and constant distribution for κ respectively and c D , c R stands for the deterministic and random distribution for c respectively. Subplot (a), (b), (c), (d) stand for eccentricity 0, h 12 , h 6 , h 3 , respectively; Figure S3. PDFs of the rotation of the settlement in rad for non porous medium. κL, κC stands for the linear and constant distribution for κ respectively and c D , c R stands for the deterministic and random distribution for c respectively. Subplot (a), (b), (c) stand for eccentricity h 12 , h 6 , h 3 , respectively; Figure S4. PDFs of the maximum settlement stresses in kPa for porous analyses with deterministic shape functions for the stochastic material variables. κL, κC stands for the linear and constant distribution for κ respectively while c D , c R stands for the deterministic and random distribution for c respectively and k R , k D stands for the random and deterministic distribution for k. Subplot (a), (b), (c), (d) stand for eccentricity 0, h 12 , h 6 , h 3 respectively; Figure S5. PDFs of the minimum displacement in m for porous analyses with deterministic shape functions for the stochastic material variables. κL, κC stands for the linear and constant distribution for κ respectively while c D , c R stands for the deterministic and random distribution for c respectively and k R , k D stands for the random and deterministic distribution for k. Subplot (a), (b), (c), (d) stand for eccentricity 0, h 12 , h 6 , h 3 respectively; Figure S6. PDFs of the rotation of the settlement in rad for porous analyses with deterministic shape functions for the stochastic material variables. κL, κC stands for the linear and constant distribution for κ respectively while c D , c R stands for the deterministic and random distribution for c respectively and k R , k D stands for the random and deterministic distribution for k. Subplot (a), (b), (c) stand for eccentricity h 12 , h 6 , h 3 respectively; Figure S7. PDFs of the maximum settlement stresses in kPa for porous analyses with random field representation for the stochastic material variables. κRF 2 , κRF 4 , κRF 8 Table S1. Monte Carlo results for the stresses (kPa), the strains (‰) and the probability of the first Gauss point failure for non porous medium. µ ν , σ δ , µ stand for the mean value, the coefficient of variation, the minimum value of the output random vector. p vol , q dev , e vol , e dev denote the volumetric stress at failure, the deviatoric Von Mises stress at failure, the volumetric strain at failure and the deviatoric scalar value for strains at failure respectively; Table S2. Monte Carlo results for the stresses (kPa), the strains (‰) and the probability of the first Gauss point failure for porous analyses with deterministic shape functions for the stochastic material variables. µ ν , σ δ , µ stand for the mean value, the coefficient of variation, the minimum value of the output random vector. p vol , q dev , e vol , e dev denote the volumetric stress at failure, the deviatoric Von Mises stress at failure, the volumetric strain at failure and the deviatoric scalar value for strains at failure respectively; Table S2. Monte Carlo results for the stresses (kPa), the strains (‰) and the probability of the first Gauss point failure for porous analyses with deterministic shape functions for the stochastic material variables.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
The following symbols are used in this manuscript:  The following tables and figures are placed in this section Table A1. Monte Carlo results for the maximum stress, (kPa) the normal force, (kN) the maximum and minimum displacements, (m) and the rotation of the foundation (rad) for non-porous medium.
µ ν , σ δ , M, µ stand for the mean value, the coefficient of variation, the maximum and the minimum value of the output random vector. σ max , N, u max , u min , φ denote the maximum footing settlement stress, the foundation axial force at failure, the maximum failure displacement, the minimum failure displacement and the corresponding rotation of the shallow foundation, respectively.             Figure A2. PDFs of the maximum displacement in m for non-porous medium. κ L , κ C stands for the linear and constant distribution for κ, respectively, and c D , c R stands for the deterministic and random distribution for c, respectively. Subplot (a-d) stand for eccentricity 0, h 12 , h 6 , h 3 , respectively. Table A2. Monte Carlo results for the maximum stress kPa the normal force kN the maximum and minimum displacements m and the rotation of the foundation rad for porous analyses with deterministic shape functions for the stochastic material variables. µ ν , σ δ , M, µ stand for the mean value, the coefficient of variation, the maximum and the minimum value of the output random vector.
σ max , N, u max , u min , φ denote the maximum footing settlement stress, the foundation axial force at failure, the maximum failure displacement, the minimum failure displacement and the corresponding rotation of the shallow foundation respectively.                       Figure A4. PDFs of the maximum displacement in m for porous analyses with deterministic shape functions for the stochastic material variables. κ L , κ C stands for the linear and constant distribution for κ, respectively, while c D , c R stands for the deterministic and random distribution for c, respectively, and k R , k D stands for the random and deterministic distribution for k. Subplot (a-d) stand for eccentricity 0, h 12 , h 6 , h 3 , respectively. Table A3. Monte Carlo results for the maximum stress kPa the normal force kN, the maximum and minimum displacements m and the rotation of the foundation rad for porous analyses with random field representation for the stochastic material variables. µ ν , σ δ , M, µ stand for the mean value, the coefficient of variation, the maximum and the minimum value of the output random vector. σ max , N, u max , u min , φ denote the maximum footing settlement stress, the foundation axial force at failure, the maximum failure displacement, the minimum failure displacement and the corresponding rotation of the shallow foundation, respectively.           Figure A7.