Probabilistic Failure Estimation of an Oblique Loaded Footing Settlement on Cohesive Geomaterials with a Modified Cam Clay Material Yield Function

In this work, a quantitative uncertainty estimation of the random distribution of the soil material properties to the probability density functions of the failure load and failure displacements of a shallow foundation loaded with an oblique load is portrayed. A modified Cam Clay yield constitutive model is adopted with a stochastic finite element model. The random distribution of the reload path inclination κ, the critical state line inclination c of the soil and the permeability k of the Darcian water flow relation, has been assessed with Monte Carlo simulations accelerated by using Latin hypercube sampling. It is proven that both failure load and failure displacements follow Gaussian normal distribution despite the excessive non-linear behaviour of the soil. In addition, as the obliquity increases the mean value of failure load and the failure displacement always increases. The uncertainty of the output failure stress with the increase of the obliquity of the load remains the same. The failure spline of clays can be calculated within an acceptable accuracy with the proposed numerical scheme in every possible geometry and load conditions, considering the obliquity of the load in conjunction with non-linear constitutive relations.


Introduction
The ultimate load of footing settlements and the subsequent displacements and material states is an area of significant concern in the field of geomechanics. The bearing capacity of soils, originally addressed by [1], has been widely researched and the literature consists of a number of articles published, such as [2][3][4][5][6]. The aforementioned publications provide a consideration of analyses with Mohr-Coulomb yield criterion and linear elastic material. In the vast majority of the cases 1D or 2D soil domains have been analyzed along with homogeneous or layered soils: [7][8][9]. The literature was subsequently applied to the foundation design regulations by introducing the shape influence variables S f and the friction variables N f . The friction variables N f are three respecting for the three addends contributing to the total bearing capacity. The influence of possible vertical load in the lateral of the foundation is controlled by the variable N q , the cohesion of the soil is controlled by the variable N c and the settlement dimensions alongside with the total weight of the soil is controlled by the variable N γ . Similarly, the shape variables S q , S c , S γ are incorporated [10][11][12].
The uncertainty estimation of the limit load of shallow foundations of saturated porous continuum with respecting to the input uncertainty, such as the Young Modulus or the permeability, has been assessed with the implementation of the stochastic finite element method, where the input variable and its stochastic spatial interpretation can be given by using two alternative methods. The first method considers the properties of nodal points following a set of random variables and the shape functions are adopted for the interpolation of the material spatial distribution as deterministic functions [4,[13][14][15][16][17]. Subsequently, random field processes may be implemented, like the spectral representation or the Karhunen-Loeve series expansion or the spatial average method [18][19][20][21][22][23][24][25]. The sampling method can be either non-biased through pseudorandom relations or an importance sampling method could be adopted, such as Latin hypercube Sampling (LHS) [26,27]. Then, the standard Monte Carlo simulation can be performed.
The scientific publications offer in the vast majority analyses that incorporate 1D-2D elastic halfspace theory with material yield function the Mohr-Coulomb criterion. In the recent years, scientific publications with 3D analyses that implement the Mohr-Coulomb yield criterion have been incorporated [28][29][30]. In the present paper, a numerical simulation scheme with a modified Cam Clay material yield model proposed in [31] is implemented, which is an accurate and valid quantitatively material for cohesive soil simulation [32]. This material yield surface with the incorporation of a finite element model can estimate the real load and displacement field in every possible 3D loading conditions. The computational cost needed by the crude Monte Carlo simulation method can be decreased with efficient computational schemes as proposed in [33,34]. For the calculation of failure load an improved version of a recurrence relation algorithm proposed by the authors of [35], which provides reliable results with a small number of trials and with one initial guess trial is portrayed. The algorithm is theoretically defined, proven and compared with the bisection method. The aim of the present paper, which is an extension of previous work in [35], is to quantify numerically, in the case of shallow foundations, the variability of the failure load and failure displacement alongside with the uncertainty quantification of the failure mechanism in relation to the variability of the input parameters like the spatial distribution of the material variables and the obliquity of the load with respect to the horizontal direction.
Three material variables are assumed to be stochastic: the compressibility factor κ, the permeability factor k and the critical state line inclination c of the soil. The compressibility factor κ can be determined by the standard isotropic compression experiment. κ is the inclination of the reload path in terms of void ratio and the natural logarithmic scale of the volumetric component of the stress tensor. The critical state line inclination c is directly associated with the friction angle and the cohesion. Consequently, c can be obtained with triaxial experiments and by implementing the Mohr-Coulomb yield criterion in the 3D principal stress plane the friction angle φ 0 and the cohesion. The permeability k can be estimated with several methods of providing a hydraulic gradient and measuring the velocity of the flow like pumping tests, falling head permeability tests and estimations through the cone penetration test (CPT). The randomness of the three material parameters can be justified through the literature [36][37][38]. For κ the randomness is evident since the reload path may vary significantly due to the structure of the clay microtiles. The randomness of c, which is a function of friction angle (φ 0 ), is also seen and may differ notably even in adjacent positions. Finally, since the water flow is not constant even with constant hydraulic gradient throughout the soil domain, k may notably diverge from its mean value.
The spatial distributions of the material parameters with respect to the depth of the soil domain are the linear, the constant variation and the random field process, formulated from the Karhunen-Loeve series expansion adopting an exponential autocorrelation function. In the constant and linear variation with respect to depth, the truncated normal random variable [39,40] is assumed at the nodal points. The LHS importance sampling method is implemented in order to obtain the input random vectors. The eccentricities of the shallow foundation in X and Y directions and correlation lengths are parametrized and compared to the analogous solid problem in which the pore pressure is neglected.

Dynamic Soil-Pore-Fluid Interaction: The System of Equations and Its Numerical Solution
The generalized porous media behaviour can be computed with a system of partial differential equations called the Biot problem. When low frequency loads and static forces are considered, the Biot problem is deducted to a more simple system which has an alleviated problem complexity and computational cost. The u-p formulation of Biot system of equations, which consists of the soil-fluid energy balance with the Darcian flow, the stress-strain constitutive relation and the boundary conditions is numerical problem which is more stable than the corresponding Biot system of equations. In the present work the u-p formulation is considered since static forces are implemented to the clay soil domain. It should be noted that for a variety of natural clays the u-p formulation is suitable in practically all the possible frequency excitations as proved in [41].
The finite element discretization of the u-p formulation takes the form following [42,43]: The augmented mass matrix M, stiffness matrix K and damping matrix C are as follows: M S is the standard mass matrix of the solid skeleton and denoting with ρ d the density of the soil and with N u the shape functions of the displacement field the formulation of the mass matrix is given by: K S is the standard stiffness matrix of the solid skeleton. B, E is the deformation and elasticity matrices respectively. The standard stiffness matrix takes the form: C S is the standard mass matrix of the solid skeleton which in the present work is a Rayleigh damping matrix. Furthermore, the load vector and the unknown variables vector are composed as: The augmented matrices consist of the following three components. The saturation matrix S = V N p 1 Q N p dv where N P are the shape functions for pore pressure and Q is a function of the bulk moduli of fluid and soil skeleton. The permeability matrix H = V ( N p ) T k N p dv where k is the matrix of permeability. The coupling matrix Q c = V B T mN p dv where m is the unity matrix. Finally, the loading vector divided by the total mixture density denoted with b, results to an equivalent force vector f S = V (N p ) T T (kb)dv. For the aforementioned numerical simulation algorithm numerical integration schemes such as the Newmark method can be used to obtain the solution.

Plastic Yield Envelope and Bond Strength Envelope Mathematical Representation
The material yield function used in this paper is a modified Cam Clay type model. In this section the effective stresses in the solid skeleton are considered. The model is defined by two surfaces, the plastic yield envelope (PYE) which defines the elastic region and the bond strength envelope (BSE) which defines the acceptable positions for PYE [31,32,[44][45][46]. BSE is influenced by the structure of the cohesive soil microtiles. If a stress point lays in the BSE boundary, the structure degradation rate of the clayey soil is maximized. Both envelopes are ellipses and are depicted in Figure 1.
The general representation of an envelope is the following: In Equation (6) 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 ξ in 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 generalized 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 to be 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.
The aforementioned material constitutive model is reliable for the majority of natural cohesive soils, of which the friction angle varies between 17 • and 30 • . The representation of the response of a clayey soil point is accurate in all possible loadings. Moreover, the equations of the criterion are in closed form which contributes to a stable numerical system. In addition, when large values of OCR of a cohesive soil is simulated this yield function can be easily reformed to take into account possible tensile stresses. The stresses and the strains algebraic transformations are performed 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 (10) The variable e is used for denoting deviatoric strain measure whilst q is used for denoting the deviatoric stress measure (Von Mises stress).

The Karhunen-Loeve Series and the Truncated Normal Variables
The material uncertainty of the random parameters can be interpreted either by considering that the nodal points following random variables with deterministic shape functions, or by using the Karhunen-Loeve series expansion for the estimation of a realization of a random field [47][48][49][50][51][52][53]. In this paper, both methods are used in order to evaluate the performance of each simulation procedure.
The method of the Karhunen-Loeve series expansion applied in the present work, is formed for the calculation of H 1 (x, ω), a random field of mean µ(x) based on the exponential autocovariance function. Denoting with b the correlation length the autocovariance function is stated as: A realization H 1 of the field considering M e number of eigenfunctions φ i with corresponding eigenvalues λ i can be computed as: where ξ i is a set of random variables of zero mean and covariance function COV(ξ i , ξ j ) = δ ij and T p is the symmetrization factor of the subspace of the process. If a stochastic process is Gaussian, as adopted in this paper, the ξ i functions are a set of standard normal random variables. This series expansion is the most widely used due to its high accuracy. If the Fredholm eigenvalue problem cannot have an analytical solution, and this takes place when the autocovariance function is more complicated, numerical methods can be introduced [19,50].
In [47] a method for interpreting the spatial randomness of a material parameter has been proposed, in which the random function f is formed with the implementation of shape functions N i . Denoting with f i the nodal point values, which are random variables following a probability density function (PDF) and N 0 the total number of shape functions the random function f takes the form: In the present work, linear shape functions are adopted and the f i follow the truncated normal distribution ( [39,40,54,55]), whose PDF is described by the equation In Equation (14), A, B, X 0 denote the normalized coordinates of the subspace limits and x, respectively while φ(X 0 ) and Φ(X 0 ) denote the standard normal probability and cumulative distribution function for X 0 . Additionally, σ d and the mean value of the normalization corresponds to the PDF of the random variable before truncation.

Latin Hypercube Sampling
The random vector samples can be selected by implementing pseudorandom unbiased algorithms to generate a large random vector, or by variance reduction methods like importance sampling [56] and Latin hypercube sampling (LHS) ( [26,27]). The choice of the random vector through LHS the method reduces notably the computational cost in order to estimate the statistical moments of the output random variables.
The Latin hypercube procedure is as follows: consider X as a random vector (x 1 ,x 2 ,. . . x n ). For each random variable x i the following are made: • The interval [0,1] of the cumulative distribution function (CDF) is subdivided into N s equal subspaces • A random number is chosen and through the inverse CDF a sample x i is acquired • All the x i vectors are permuted in a random way and thus the vector realization X is composed The aforementioned procedure ensures that at each possible row and each possible column in the (n x n) Euclidean space one sample is exactly extracted. This attribute is portrayed in Figure 2. The most profound advantage of this scheme is that a smaller amount of values in comparison with the crude Monte Carlo simulation are required to integrate the PDF of the input uncertainty and subsequently to approximate the variability of the output. Moreover, the subintervals in each dimension are of different sizes in order to take into consideration possible asymmetries of the PDF of the input. This perspective can be used to cross-correlate variables, i.e., when the correlation matrix is not diagonal, thus it is of general application in the uncertainty quantification in engineering systems. In the present work the variables are not cross-correlated and the LHS sampling is implemented to acquire samples for the variables in the 3-dimensional space of the compressibility factor κ, the critical state line inclination c and the permeability k.
The material random variables following the PDF h 1 of Section 4.1 and the random field realizations calculated with the Equation (12) effect the finite element numerical scheme of Equation (1). The corresponding matrices C, K and F alter due to the randomness of the compressibility factor κ, as well as of the critical state line inclination c and the permeability k. The selection of the samples follows the importance sampling technique of the LHS method.

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 of [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 3. The objective of this algorithm is to estimate the load factor λ * which causes failure of a body at exactly the time T. An initial trial of λ 1 is considered leading to an initial time of failure t 1 is obtained. Then, for each new estimation λ n+1 if the load factor λ n causes failure to the body is calculated by the equation If the load factor λ n is not causing failure, then the maximum no failure factor λ max−no− f ailure is acquired from all previous estimations used for the calculation of λ n+1 , with the implementation of the following equation: In practice this recurrence relation usually converges "by the failure region" which means that only Equation (15) is implemented. The difference between λ n+1 and λ n is given by: In Equation (17) it is evident that as n → ∞, t n → T consequently, the left side of the equation tends to 0, and consequently the algorithm converges to the desired load factor λ * .
In Table 1 an investigation of the proposed algorithm is depicted. Comparing with the typical bisection method, where one should estimate initial values of failure and safety values λ 1, f ail and λ 1,no− f ailure and then compute the new load factor by the bisection of maximum safety factor and minimum failure factor, this algorithm gives less or an equal number of trials for the same initial failure guess and convergence tolerance, as is portrayed in Table 1. In addition, it needs only one initial guess, which makes the divergence of the solution less probable. In conclusion, as proven by the numerical tests presented in Table 1 and in Section 6, the divergence in failure load and failure displacement between the two algorithms is less than 1%. The absolute percentage difference is computed considering as exact solution of the bisection algorithm with an initial value of safety stress 1000 kPa. The computational performance for the 100 deterministic analyses of a Monte Carlo simulation is also presented in Table 1. The proposed algorithm can save up to almost 55% of the time demanded for a Monte Carlo analysis to be performed, indicating a considerable advantage of the aforementioned recurrence relation.

Description of the Problem
The proposed numerical simulation scheme is implemented in porous problems as portrayed in Figure 4 and are defined by the set of Equations (1). The monitored output variables are the total force of the footing settlement and its maximum and minimum displacements which are the mean value of displacements calculated at points A-D of Figure 4. The loading conditions consist of the oblique nodal values q 1 − q 4 at points B, C, A and D respectively and incorporate the equivalent forces of an oblique shallow foundation load having angle of obliquity denoted as θ q , with respect to the horizontal direction and is uniformly distributed over the footing settlement surface. The modelling of the foundation consists solely of its equivalent forces along its surface ABCD, which is measured (1 × 1 m 2 ). The finite element mesh is constructed with eight node hexahedral finite elements with linear shape functions for u and p, which gives an acceptable numerical reliability [57,58]. The soil dimensions in X, Y, Z directions are respectively l x = 5 m, l y = 5 m, l z = 4 m. The geostatic stresses are directly imported as initial conditions, with the relations σ v = γz, σ x = σ y = 100 kPa and are associated with stress point L of Figure 1. The simulation duration in all cases is one day and quasi static conditions are attained. A time step of dt = 0.001 d is adopted. The other deterministic parameters of the soil are given in Table 2. In Table 2, λ denotes the inclination of isotropic compression line for the respective virgin normally consolidated clay and ν 0 denotes the initial specific volume of the soil, which is assumed proportional to κ. The boundary conditions are: u x (z = h) = u y (z = h) = u z (z = h) = 0 and the lateral boundary surfaces are constraint-free. The input material uncertainty consist of the material variables, the compressibility factor κ, the critical state line inclination c and the permeability factor k.
The compressibility factor κ, is considered as constant (κ C ) over depth, or linear (κ L ). In the κ L occasion, κ 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 chosen in order for the mean stiffness of the soil to correspond to a shear velocity of 200 m s . In Figure 5 the spatial distribution of κ L is portrayed. It should be emphasized that that bulk and the shear moduli are proportional since the Poisson ratio is assumed as constant. As a consequence, κ is directly related to the shear velocity. When κ has a constant distribution over depth, the mean value of κ is κ µ = 0.004074 and the CoV is 0.25.  For the critical state inclination c the constant variation over depth is presumed. Two possible considerations for the absolute value are implemented. If a random variable case is assumed c R , the friction angle φ 0 follows the truncated normal distribution with the mean value is µ φ = 23 • and the standard deviation is σ φ = 2 • which gives acceptable values for φ 0 considering for natural clays ( [31]). A random vector of φ 0 is obtained through Latin hypercube sampling for the standard normal distribution. The vector φ 0 components are transformed into the truncated normal PDF h 1 . Finally, c is computed from c = 2 3 6sin(φ 0 ) 3−sin(φ 0 ) . In an alternative case, for the deterministic value c D , c = 0.7336 for friction angle µ φ = 23 • .
The permeability k, is considered constant over depth. The absolute value may be calculated by two possible cases. If a random variable approach 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 investigated. The solid analyses, where the water flow is not taken into account and the porous analyses, where the pore pressure is considered. The solid analyses performed, denoted with (S), are presented in Table 3, using linear (L) or constant (C) distribution for κ and deterministic (D) or random variable (R) cases for c. The corresponding porous analyses performed abbreviated are depicted 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. Deterministic (D), random variable case (R) and random field distribution (RF) for k.
In the stochastic processes the mean values are chosen as follows: κ mean = 0.008686, c mean = 0.7336 and k mean = 10 −8 m 3 s Mgr following [41,43,59]. The standard deviations chosen are: σ κ = 0.25κ mean , σ φ = 2 • and σ k = 0.25k mean . The exponential autocorrelation function described in section Section 4.1 is used in all stochastic processes. The correlation lengths are chosen b = 2 m (k RF ), b = 4 m (k RF ) and b = 8 m (k RF ). The spatial distributions for κ, κ L and κ C , as well as the random variable distributions for all material variables are a random variable case analysis. For c a constant deterministic analysis is assumed. 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 Section 4.1 with the use of the autocovariance function given by Equation (11).
The simulations are static, and the number of Fredholm eigenfunctions taken into account 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 implemented with 100 samples, using the Latin hypercube sampling method, which led to convergence for the mean value and standard deviation of the monitored displacements as depicted in Figure 6 where 1000 samples were taken for a randomly selected Monte Carlo simulation and compared to the corresponding statistics for 100 samples. It should be mentioned that the cross-correlation of all material variables is not taken into account and subsequently the correlation matrix is diagonal.

Limit Failure Force and Displacement
The results for the statistics of failure force and failure displacement are presented in Tables A1-A3 and in Figures A1-A9 in Appendix A. In these tables the mean values (µ ν ) of the total settlement force N, the vertical displacement u y and the horizontal displacement u x at failure are portrayed with the coefficient of variation (CoV) (σ δ ) as well as the maximum (M) and minimum µ values acquired from the Monte Carlo simulation. The PDFs of the simulations are portrayed in the presented figures. They have been obtained by implementing the Kolomogorov-Smirnov test for proving their Gaussian normal distribution and then the corresponding histograms are smoothened for drawing the probability density function.
As portrayed in Table A1, when the water flow is not present, larger mean failure displacements and lower CoV are obtained when κ L is assumed compared to κ C . In addition, a smaller mean total foundation load is obtained when κ L case is adopted whilst the corresponding uncertainty is maximized for the κ c − c R case. In solid analyses the maximum CoV of the output for failure footing force is 40% the input variability while the largest CoV for failure displacement in both directions is the same as the input uncertainty (see κ C − c R for all monitored variables). This applies for all values of obliquity angle θ q . The mean value of N for θ q = 60 degrees is 1,9 times than the value for θ q = 0. The mean values for the u y increase with the increase of θ q while for u x the critical angle for the maximization of its mean value is θ q = 30. Thus, the critical spatial distribution of κ for the CoV of the output in all monitored variables is the κ C case and for c, the random variable case. The PDFs of solid analyses are portrayed in Figures A1-A3 where it can be concluded that the increase of obliquity leads to an increase of N and u y while for u x for obliquity angles of more than 30 degrees, a significant decrease is observed. This behaviour in the displacements is explained by the fact that in κ L distribution the upper compressible layers of the soil, have smaller CoV of κ causing smaller variability for displacements and strains. The assumption of constant variation over depth of the soil mass tends to be more homogeneous and stiff and, consequently the Gauss points have greater stiffness and larger failure footing forces occur.
The output uncertainty of the foundation total load in porous analyses is influenced to the same magnitude compared to solid problems with the alteration of obliquity. The variability of the output failure displacements increases with the increase of the obliquity until the angle of 45 degrees and then decreases, unlike the corresponding solid analyses in which the obliquity has practically no influence on the CoV of the output. As depicted in Table A2, the largest CoV of the footing total load, which is found at P−κ C − c R − k d , θ q = 60 • , is 40% of the variability of the input, while for both horizontal and vertical displacements the maximum CoV is 12% greater than the uncertainty of the input and is located at Therefore, when taking into account the pore pressure in the soil domain, a decrease of the variability of the failure footing force takes place whilst in failure displacements when there is constant variation over depth for κ, a significant uncertainty enlargement occurs, as indicated in Figures A4-A6. This way of behaving can be attributed to the fact that when a constant spatial variation of a material variable is adopted, all Gauss points are associated with greater uncertainty leading to larger CoV for strains and displacements. Taking into account that the bulk modulus is related through Equation (9), in porous problems it is generally smaller than the corresponding solid problems, it is concluded that smaller mean values of failure force and smaller variability are expected, due to tensile failure of the first Gauss point. This will be proven numerically in Section 6.2.2 in Tables A4-A6. Porous analyses with the assumption that the material variables follow a stochastic process are conducted as a more general case in order to take into account the spatial randomness of the material parameters of the soil. In porous random field simulations, the greatest CoV of the shallow foundation total load is 48% the input uncertainty, while for the displacements the reduction of the variability is by 40% and 56% for the vertical and horizontal displacement, respectively. As depicted in Table A3, the critical correlation length for maximizing the output uncertainty is b = 2 m whilst for minimizing the mean foundation force is b = 4 m for all obliquities examined. When b = 2 m, the vertical mean failure displacements are the largest whilst for mean horizontal displacements for θ q = 0 the critical value for b = 2 m and when θ q > 0 the corresponding critical failure displacements occur when b = 8 m. The mean values of failure total foundation force in porous random field analyses are slightly smaller, with divergence on the vicinity of 10% compared to porous analyses with deterministic shape functions for the material parameters κ, c, k while the mean values for the failure displacements in both horizontal and vertical directions are larger in the porous analyses with random field representation for all the material variables. In the vertical displacements the enlargement compared to the corresponding porous analyses with deterministic shape functions for κ, c, k is in the order of magnitude of 20% whilst in the horizontal displacements it is in the vicinity of 30%. The unfavourable situation takes place when the foundation force is small and, consequently, the critical spatial distribution for the mean value of N is the random Karhunen-Loeve expansion for all material input variables and b = 4 m. The probability density functions of P-κ RF -c RF -k RF are portrayed in Figures A7-A9.
The results acquired give qualitative statements and quantitative proposals on the influence of the input uncertainty of each material variable in porous failure phenomena taking into account the obliquity of the shallow foundation. The compressibility factor κ affects the statistical moments of all monitored variables. This effect in CoV of the output is more profound when κ is uniformly distributed along the depth. This holds in both porous and non-porous problems which can be attributed to the fact that κ is directly associated with the bulk modulus and consequently this reflects on the strains the displacements and the stiffness of the soil domain leading to larger alteration of the footing force.
The permeability k has alleviated the effect of failure shallow foundation forces and failure displacements. The spatial variability of k appears to influence to a lesser extent the CoV of the output in most Monte Carlo analyses. For porous consolidation problem with the load geometry and deterministic parameters remaining the same, the output displacements and stresses field are not influenced by the permeability since the pore pressures are fully dissipated, as verified by the numerical results.
Finally, the uncertainty of a plastic variable like the critical state line inclination c of the material model appears to have the most notable effect on the uncertainty of the shallow foundation total load and the corresponding displacements when the constant variation over depth and random variable case c R is adopted. The mean values for failure force are critical when c R case is assumed. For mean failure displacements the alteration of c has a negligible effect in all obliquities in both solid and porous problems. This can be explained by the fact that c is directly influenced by the integration point state at failure, subsequently it alters the limit load of the soil domain.
The validity of the PDFs presented has been verified in two ways. Firstly, the histograms obtained from each Monte Carlo simulation are smoothed and the aforementioned curve is estimated as a Gaussian normal PDF or a Gaussian lognormal PDF respectively. The Gaussian lognormal PDF is estimated when the CoV is large. Subsequently, the Kolomogorov and Smirnov test for the numerical validation of this assumption is performed. In all cases the null hypothesis for the 5% reliability interval is accepted so in all cases either the normal PDF or the lognormal PDF is thpropriate for the output monitored variables. In the present work since the largest CoV of the output obtained provide a small probability of negative values for the monitored output variables, the normal probability density functions are assumed.

Limit Stresses-Strains and Failure Spline
The results for the limit stress-strain statistics and failure splines are portrayed in Tables A4-A6. In these tables the following statistics are given: mean values, CoV and minimum values for the volumetric stresses denoted as p vol , and deviatoric Von Mises stresses denoted as q dev , as explained in Section 3.1. Moreover, the mean value with the CoV of the volumetric strains denoted as e vol and the deviatoric strains denoted as e dev are also presented in Tables A4-A6, The probability of first Gauss point failure for all simulations is also depicted. At this point, it is declared that the volumetric strain at failure is tensile.
As depicted in analyses of non-porous medium of Table A4, larger uncertainty of the output and reduced minimum value of failure volumetric and deviatoric stresses are acquired when c R case is assumed. This assumption for critical state line inclination give the greatest variability of the strains at failure in all obliquities. The mean value for the volumetric component is in the vicinity of 15‰ for the linear variation over depth for κ and in the vicinity of 7‰ for the constant variation over depth for κ. The respective values for the deviatoric component are 18‰ and 10‰. From the output results of the Monte Carlo simulation it can be concluded that the percentage of plastic deviatoric strains are larger when θ q = 0.30, thus the distortional failure is critical. For θ q = 45.60 the volumetric failure is evident. Finally, when the obliquity is not present, which means when pure horizontal loading is applied, (3.21, 2.21, 3.79) is the critical integration point, whilst when θ q > 0, (1.79, 2.21, 3.79) is the critical Gauss point.
In porous simulations with deterministic spatial distributions for all the material parameters larger stress variability in both volumetric and deviatoric components is obtained when the κ L − c R case is assumed, as portrayed in Table A5. The alteration of the permeability significantly affects the output uncertainty when θ q = 60 degrees in P−κ C − c R − k D and P−κ C − c R − k R analyses. The numerical Monte Carlo simulations indicate that the largest output uncertainty for volumetric strains is found at θ q = 60, P−κ C − c R − k R while for deviatoric stresses is located at P−κ C − c R − k D for obliquity angle of 60 degrees. The mean value of volumetric strains varies between 5 and 8‰ and the mean values of deviatoric strains is in order of magnitude of 10‰, considering constant variation over depth for the compressibility factor. For κ L case the mean values are notably larger. From the output results attained it is evident that for small obliquities the volumetric failure is critical, while for larger angles, with the exception of θ q = 60 and P−κ C − c R − k D case the distortional failure is evident. Finally, the most probable failure Gauss point is the (3.21, 2.21, 3.79) for θ q = 0, (3.79, 2.21, 3.79) for θ q = 30 and (1.79, 2.21, 3.79) for θ q = 45, 60. As the obliquity increases more, integration points share a larger probability of failure. The smallest probability of the most profound integration point is over 80% indicating the importance of this Gauss point in the estimation of the failure spline.
In porous analyses with random process for κ, c, k, greater output uncertainty for stresses is obtained for b = 2 m, in all cases as portrayed in Table A6. The volumetric strains in most obliquities have critical correlation length b = 4 m and similar conclusions can be made for the deviatoric strains. Taking into account the numerical results obtained by the analyses simulated when b = 2 and θ q = 30 degrees the probability is divided to larger amount of integration points, consequently the variability of the failure mechanism is maximum. When the obliquity angle in degrees is θ q = 0.30 the volumetric failure occurs whilst if θ q = 45.60 the deviatoric failure takes place. Finally, when θ q = 0 the critical integration point is (3.21, 2.21, 3.79), when θ q = 30 the critical Gauss point is (3.79, 2.21, 3.79), when θ q = 45, 60 the aforementioned spline point is (1.79, 2.21, 3.79).

Conclusions
In this paper, the variability estimation of the failure of shallow foundations on cohesive soils taking into account the pore pressure soil interaction with the use of the stochastic finite element method is depicted. The aim of this article is to provide a numerical tool for acquiring accurate quantitative results on the failure footing force and corresponding displacements in relation to the input variability of soil material parameters. The aforementioned numerical simulation model is valid for general assumptions for the geometry, the loading and the material distribution of the soil mass. In this context, a detailed finite element simulation, alongside an accurate and reliable material constitutive yield function are used.
The numerical results of the simulations portray that the monitored output variables of ultimate foundation load and corresponding displacements follow a Gaussian random distribution despite material non-linearity of high magnitude that is profound in failure phenomena. The randomness of material poroelasticity is important in the uncertainty of the output variables especially when it has a constant variation with respect to the depth of the soil domain. Analogous considerations can be made for the critical state line inclination c. When the constant distribution for κ is adopted, the CoV of the output maximum displacement is increased in relation to the input variability by 12%. The porous variable of permeability appears to have an alleviated influence to the uncertainty of the monitored variables. In porous problems the failure load when the obliquity is not 0, the mean failure foundation force is smaller in comparison with the corresponding non-porous medium.
The random fields for κ, c, k give larger mean failure volumetric and deviatoric strains for all obliquity angles examined. In porous stochastic process simulations when the load angle with respect to the horizontal direction is small, the volumetric failure is critical whilst when the obliquity angle is larger the distortional failure is critical. If deterministic shape function over depth is assumed for κ, c, k, similar conclusions can be made. In conclusion, when an a horizontal force is applied the

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

Nomenclature
The following symbols are used in this manuscript:

Appendix A
The following tables and figures are placed in this section.                                  Table A4. Monte Carlo results for the stresses (kPa), the strains (‰) and the probability of the first Gauss point failure for non-porous medium (θ q in degrees).              Table A5. Monte Carlo results for the stresses (kPa), the strains (‰) and the probability of the first Gauss point failure for porous analyses with linear and constant distribution for κ and c (θ q in degrees).