Probabilistic Bearing Capacity of a Pavement Resting on Fibre Reinforced Embankment Considering Soil Spatial Variability

This paper intends to examine the influence of spatial variability of soil properties on the probabilistic bearing capacity of a pavement located on the crest of a fibre reinforced embankment. An anisotropic random field, in combination with the finite difference method, is used to carry out the probabilistic analyses. The cohesion and internal friction angle of the soil are assumed to be lognormally distributed. The Monte Carlo simulations are carried out to obtain the mean and coefficient of variation of the pavement bearing capacity. The mean bearing capacity of the pavement is found to decrease with the increase in horizontal scale of fluctuation for a constant vertical scale of fluctuation; whereas, the coefficient of variation of the bearing capacity increases with the increase in horizontal scale of fluctuation. However, both the mean and coefficient of variation of bearing capacity of the pavement are observed to be increasing with the increase in vertical scale of fluctuation for a constant horizontal scale of fluctuation. Apart from the different scales of fluctuation, the effects of out of the plane length of the embankment and randomness in soil properties on the probabilistic bearing capacity are also investigated in the present study.


INTRODUCTION
Soil reinforced with fibres in optimum quantity is proven to be one of the most efficient and economical means of ground improvement technique over the past few decades. Fibre reinforcement not only improves the unconfined compressive strength and shear strength properties of the soil (Cai et al., 2006;Consoli et al., 2010), it also increases the tensile strength (Tang et al., 2016;Cristelo et al., 2017). Both the natural and synthetic fibres are used as reinforcement to enhance the mechanical properties of the soil. Because of cost-effectiveness, easy availability, and eco-friendly nature, the use of natural fibres (such as jute, sisal, coir, etc.) has gained popularity over the period. Several studies are there on the soil reinforced with natural fibre (Prabakar and Sridhar, 2002;Ghosh et al., 2005;Babu and Vasudevan, 2008;Chaple and Dhatrak, 2013;Singh and Bagra, 2013). The major drawback of using natural fibre is its biodegradation caused by microorganisms in soil which may reduce the long term applicability of the reinforcement (Tang et al., 2016). Thus it gives rise to the use of synthetic fibres as soil reinforcement. Many researchers have studied the effectiveness of using synthetic fibre-reinforced soil (Michalwoski and Cermak, 2002;Kumar et al., 2007;Tang et al., 2007;Consoli et al., 2010;Chen et al., 2015;Correia et al., 2015;Ates, 2016;Bouaricha et al., 2017;Cristelo et al., 2017;Sharma and Kumar, 2017). The reinforcing effects of fibres on the behavior of soils have been investigated numerically by several authors (Babu et al., 2007;Toh et al., 2017;Arora and Kumar, 2019;Sharma and Kumar, 2019;Wang et al., 2019).
All of the studies, as mentioned above, are deterministic in which the soil is assumed to be a single homogeneous layer or a stratified medium with uniform soil properties. It is well known that soil is very much heterogeneous and random because of its different geological formation processes and mineralogical constituents. Thus considering soil as a homogeneous medium may lead to the unreliable design of civil engineering structures. Many researchers conducted the probabilistic studies assessing the randomness in the soil as well as spatial variability (Fenton and Griffiths, 2001;Griffiths et al., 2002;Fenton and Griffiths, 2005;Luo et al., 2011;Cassidy et al., 2012;Li et al., 2014;Jha, 2016;Zaskórski et al., 2017). There is a huge possibility that the mixing of fibres with soil may not be uniform, which may lead to the randomness in fibre reinforced soil. Very few probabilistic and reliability based studies are available on fibre reinforced soil (Ranjan et al., 1996;Moghal et al., 2016;Diab et al., 2017;Moghal et al., 2017;. Johari and Kalantari (2016) carried out a probabilistic analysis of slope stability of embankment reinforced with discrete polypropylene fibre. But spatial variability of soil properties was not considered in their study. However, it is essential to consider the spatial variability as over a certain length of soil domain the strength properties may change. In past, many literature (Griffiths and Fenton, 2004;Suchomel and Masin, 2009;Kasama and Whittle, 2015) considered the soil spatial variability in slope stability problems. Recently, Luo and Bathurst (2018); Halder and Chakraborty (2019) have carried out the probabilistic analyses on geogrid reinforced embankment considering soil spatial variability to investigate the loadsettlement behavior.
The main objective of this paper is to investigate the effect of soil spatial variability on the performance of a pavement located on the top of the fibre-reinforced embankment. Different horizontal and vertical scales of fluctuation are chosen to study the influence of soil spatial variability. It is to be noted that all the probabilistic studies stated above were based on plane strain condition. To the best of authors' knowledge, there is no probabilistic study available on the fibre-reinforced embankment considering spatially variable soil parameters. Note that, for spatially variable soil it is essential to consider the length of the model in out of plane direction as it significantly affects the mean and coefficient of variation of load-carrying capacity of the geotechnical structures (Kawa and Pula, 2019). Hence, the influence of different out of the plane lengths of the embankment on the bearing capacity of the pavement is also explored based on both deterministic and probabilistic studies. Here only the embankment soil is assumed to be fibre-reinforced. In contrast to that, the foundation soil is kept as unreinforced. The randomness and spatial variability of the soil properties are considered for both embankment and foundation soil, and their effects are included by associating the finite difference mesh and random field. Three-dimensional finite-difference software (FLAC 3D ) is used to carry out the numerical analyses. The probabilistic results are obtained using the Monte Carlo simulation technique. Finally, failure probabilities of the bearing capacity of the pavement are computed for different horizontal and vertical scales of fluctuation. Figure 1 depicts the embankment slope geometry considered in the present numerical analyses. The embankment having a height of 5 m and an inclination of 1.5 H:1 V is filled with fibrereinforced soil consisting of improved cohesion (c) and angle of internal friction (ϕ); whereas, the foundation soil is assumed to be unreinforced with lower shear strength properties. A pavement of width 3.5 m (B 3.5 m) is located at the crest of the embankment. The present work means to estimate the probabilistic bearing capacity of the pavement for a certain settlement value, which can be expressed as mean (µ q ) and coefficient of variation (COV q ) of bearing capacity of the pavement.

NUMERICAL MODELING
Three-dimensional explicit finite difference software, FLAC 3D is engaged to generate the three-dimensional modeling of the embankment and to accomplish the numerical analyses. The domain size of the problem, in both horizontal and vertical direction, is chosen in such a way that there should not be any boundary effect. The displacement along the bottom boundary edge is fixed in both horizontal and vertical directions; whereas, the side boundaries are horizontally restricted in order to allow the vertical displacement only. Eight noded brick shaped elements are used for discretization of the problem domain. A comparatively finer mesh is generated to model the embankment soil; while, coarser mesh is chosen in order to model the foundation soil. The Mohr-Coulomb yield criterion is incorporated to simulate the behavior of embankment and foundation soil. However, it should be mentioned that embankment soil and foundation soil may differ in terms of type, composition and strength requirements. Hence, the use of only Mohr-Coulomb yield criterion for both type of soil may not represent the real scenario. It can be considered as one of the limitations of this study. The soil parameters of fibre-reinforced embankment soil and unreinforced foundation soil considered in the numerical analyses are taken as provided by Sharma and Kumar (2019) and given in Table 1.
After modeling the embankment slope geometry and allocating the soil properties node wise, the loading is simulated on the nodes representing pavement width (and along with the out of the plane length of the pavement) by implementing a very small amount of downward velocity in vertical direction. An optimized velocity of magnitude 5 × 10 -6 m/step is taken after few trials, as it is found to be less time-consuming as well as does not affect the pavement bearing pressure settlement response (Halder and Chakraborty, 2020). The numerical model is then run for numerous steps until the plastic steady state is achieved.

VALIDATION OF THE PRESENT STUDY
Before carrying out the numerical analyses, the finite difference model is validated with the available study. Since there is no probabilistic study available on fibre-rerinforced embankment considering soil spatial variability, the present result is verified with the deterministic results of Sharma and Kumar (2019). They carried out a deterministic three-dimensional numerical study on bearing capacity of ring and circular foundation resting on twolayered sand using finite element method in which top layer sand is reinforced with fibre, and the underlying layer of sand is kept as unreinforced. In the present work, only the circular footing and the two-layered soil system are modeled with identical geometry and soil properties. The bearing capacity and footing settlement curve is obtained and compared with that of Sharma and Kumar (2019). Figure 2 illustrates that the present bearing pressuresettlement response is comparatively lower than that obtained in literature which may be due to the difference in mesh generation and applied numerical scheme. The obtained result from present study is observed to be in the conservative side.

DETERMINISTIC ANALYSIS
Before implementing the probabilistic analyses, deterministic analyses are carried out on the numerical model where both the foundation and embankment soil properties are considered to be uniform or homogeneous. At first, the whole system is modeled using unreinforced soil properties for different out of the plane lengths (L op ). Then only the embankment soil is modeled considering reinforced soil properties. Both the unreinforced and reinforced bearing capacity-settlement response curves obtained for different L op are illustrated in Figure 3. It is quite obvious that the embankment soil reinforced with fibre is having a higher bearing capacity than  the unreinforced soil. As an example, the bearing capacity of the pavement for L op 3 m is increased from 17.5 kPa to 56.6 kPa corresponding to 60 mm of settlement value when the unreinforced embankment is replaced with fibre reinforcement. The bearing capacity of the pavement for unreinforced embankment tends to increase from 17.7 kPa for L op 0.5 m to 17.8 kPa for L op 1 m, then it decreases to 17.5 kPa for L op 3 m and further increase to 17.7 kPa for L op 5 m corresponding to the same settlement value as stated earlier. In case of reinforced embankment, it decreases from 56.9 kPa for L op 0.5 m to 56.6 kPa for L op 2 m. Beyond L op 2 m, it remains almost unaltered. Although the differences in bearing capacity for different L op are very negligible, it signifies the importance of carrying out the probabilistic analysis to investigate the influence of different L op on the probabilistic bearing capacity of the pavement.

Random Field Generation
The variations of the properties of the in-situ soil can be represented by the mean value, coefficient of variation and scales of fluctuation (SOFs) (Haldar and Babu, 2008). Most of the studies (Griffiths et al., 2002;Haldar and Babu, 2008;Ahmed and Soubra, 2014) have considered the cohesion (c) as lognormally distributed random field represented by mean (µ c ) and standard deviation (σ c ). The lognormal distribution is chosen to avoid the generation of negative values of soil parameters. Due to the fundamental nature of the parameter tanϕ in the equation of Mohr-Coulomb yield criterion, it is modeled as a lognormally distributed random field rather than ϕ itself. A lognormal distribution of tanϕ ensures that the friction angle is bounded by 0 < ϕ < 90° (Griffiths et al., 2011). A lognormally distributed random field can be expressed as: whereX X (x,ỹ,z) is the spatial position at which c and ϕ are desired. G(X) is a normally distributed random field with zero mean and unit variance. The values of μ ln c , μ ln tanϕ and σ ln c , σ ln tanϕ are determined using Lognormal distribution transformation given by where μ tan ϕ and σ tan ϕ are the mean and standard deviation of tanϕ, which is lognormally distributed. The correlation function [ρ(τ)] which is also known as Markov correlation function, can be expressed as where τ x |x 2 −x 1 |, τ y ỹ 2 −ỹ 1 and τ z |z 2 −z 1 | are the absolute distance between two points. Parameters δ x , δ y and δ z are the SOFs in x, y and z directions, respectively. The correlation matrix is decomposed into the product of a lower triangular matrix (L t ) and its transpose by Cholesky decomposition, Using the lower triangular matrix, the random field can be generated which is shown by where Z j is the sequence of independent standard normal random variables.
In the present study, both the c and tanϕ for embankment soil are considered as random variables. Since the foundation soil is having very low value of cohesion (c 1 kPa), it is considered as constant; whereas, only tanϕ is chosen as the random variable.
The mean values of c and tanϕ used in the probabilistic analysis are taken as the constant property values for deterministic analysis. Typical values of coefficient of variation (COV) for tanϕ and c and the horizontal SOF (δ x /B δ y /B) as well as the vertical SOF (δ z /B) are selected from Phoon and Kulhawy (1999). The probabilistic parameters considered in the present study are listed in Table 2. It should be mentioned here that there is no cross-correlation considered between tanϕ and c.
The nodal coordinates of the finite difference mesh are taken from FLAC 3D and imported to MATLAB. In MATLAB, the Markov correlation function given in Eq. 7 is used to generate the random field. The randomly distributed c and ϕ values generated in MATLAB are again taken back to FLAC 3D and are allocated node wise (Halder and Chakraborty, 2018;Halder and Chakraborty, 2020). In this way, the random fields are generated in FLAC 3D .

Monte Carlo Simulations
For each set of statistical parameters, such as COV c , COV tanϕ and δ i (where i x, y and z), the Monte Carlo simulations are carried out in probabilistic analyses to evaluate the mean and COV of the bearing capacity (µ q and COV q ) of the pavement. It should be mentioned here that the number of realizations of the Monte Carlo simulations should be such that the stable solution of µ q and COV q are achieved. The fluctuations between two consecutive realizations of µ q and COV q should fall within a tolerable range, which is between 5% and 10% (Haldar and Babu, 2008

Results Obtained in the Probabilistic Analyses
Effects of the spatial variablity, randomness in the soil properties and different out of the plane lengths of the embankment on the mean and COV of bearing capacity of the pavement are discussed in the following sub-sections. All the probabilistic analyses carried out here are only for fibre-reinforced embankment. As the variation of soil properties of natural deposit in the horizontal direction is generally quite less compared to that of the vertical direction due to the process of deposition, an anisotropic random field is generated for the present study where the vertical SOFs are chosen to be less than the horizontal one. The horizontal and vertical SOFs are varied as per Table 2. However, during the parametric study, there may be some situations where isotropic random fields are generated and vertical SOF is greater than horizontal SOF.
Effect of the out of the Plane Length of the Embankment (L op ) on Probabilistic Bearing Capacity The effect of different L op on µ q and COV q is investigated by considering δ x /B δ y /B 2 and δ z /B 0.5 for both embankment and foundation soil, and COV c 25% for embankment soil only. The effects of two different coefficients of variation of soil friction angle (COV tanϕ ) for both embankment and foundation soil are also studied in this section. Figures 5A,B illustrate the behavior of µ q and COV q for different L op and COV tanϕ . As per Figure 5A the magnitude of µ q for a particular value of COV tanϕ , at first tends to decrease with an increase in L op from 0.5 to 1 m, then it increases and after L op 3 m the change in µ q is found to be insignificant.
As an example, µ q decreases from 52.36 kPa for L op 0.5 m to 51.95 kPa for L op 1 m and then increases to 52.41 kPa for L op 3 m for a constant value of COV tanϕ 20%. The increase in µ q beyond L op 3 m is quite less. The effect of COV tanϕ on µ q is quite prominent. It is evident that with an increase in COV tanϕ , there is an increase in variability in the angle of internal friction, which in  turn reduces the value of µ q . For a particular value of L op 3 m, µ q decreases from 54.43 to 52.41 kPa with an increase in COV tanϕ from 10 to 20%. However, in all the cases, the magnitudes of µ q are less than the deterministic bearing capacity values. The effect of L op on COV q is depicted in Figure 5B. In contrast to µ q , COV q is decreasing with the increase in L op . As an example, COV q decreases from 8.45 to 6.6% as L op increases from 0.5 to 5 m for a constant value of COV tanϕ 20%. Unlike µ q , COV q increases with the increment in the magnitude of COV tanϕ as the increase in COV tanϕ causes an increasing randomness of angle of internal friction. For example, the COV q increases from 4.25 to 7.06% with the increase in COV tanϕ from 10 to 20% for L op 3 m.

Random Field Plots Considering Spatial Variability
The random fields are generated by implementing the aforesaid theory using FLAC 3D and MATLAB. Since there is no significant difference in µ q beyond L op 3 m, all the studies considering spatial variability are done for L op 3 m only. Figures 6A-D represent the random field plots of ϕ for both embankment and foundation soil with different δ x /B δ y /B, δ z /B, COV tanϕ 20% and out of the plane length (L op ) 3 m for a particular realization. Figures 6A,B show that randomness of ϕ decreases in horizontal direction with the increase in δ x /B δ y /B value from 0.5 to 4 for the constant value of δ z /B 0.5. Similar observation can be made for variation of ϕ in vertical direction from Figures 6C,D where δ z /B is increased from 0.25 to 2 for the constant value of δ x /B δ y / B 2. Figures 6A,D  The δ x /B δ y /B, δ z /B and COV tanϕ are kept as same for both embankemt and foundation soil; whereas, COV c 25% is considered for embankment soil only. The lower value of horizontal SOF indicates that the soil friction angle field is very much erratic in horizontal directions whereas the increased value of horizontal SOF specifies the uniform nature of the soil friction angle in horizontal direction. Figure 8A depicts that the µ q decreases with the increase in δ x /B δ y /B values. In contrast to that, the COV q increases with δ x /B δ y /B as illistrated in Figure 8B. For an instance, the values of µ q decreases from 53.78 to 51.91 kPa with the increase in values of δ x /B δ y /B from 0.5 to 4, the reduction in µ q beyond δ x /B δ y /B 4 is quite negligible; whereas, the values of COV q increases from 3.73 to 10.55% with the increase in values of δ x /B δ y /B from 0.5 to 10. However, the rate of increase in COV q beyond δ x /B δ y /B 4 is comparatively less than that of before δ x / B δ y /B 4.

Effect of Vertical SOFs on the Probabilistic Bearing Capacity
The influences of different vertical SOFs (δ z /B) on the probabilistic bearing capacity of the pavement are demonstrated in Figures  9A,B for the constant values of L op 3 m, COV tanϕ 20%, COV c 25%, and δ x /B δ y /B 2. In the present study, the values of δ z /B are varied from 0.25 to 8. Unlike the horizontal SOF, the µ q is found to be increasing with the increasing magnitude of δ z /B. The obtained trend is also true for COV q . For an example, the value of µ q increases from 52.03 to 54.57 kPa and COV q increases from 5.66 to 10.8% as the δ z /B increases from 0.25 to 8. However, the rate of increase in µ q and COV q beyond δ z /B 2 is considerably less than that of before δ z /B 2.

Effect of Soil Spatial Variability on the Failure of the System
In the present work, the failure of the pavement is presented through the maximum shear stress contour profiles. Figures 10A-C demonstrate maximum shear stress contours for the embankment as well as the foundation soil system having either homogenous soil field or spatially distributed soil field with L op 3 m. The ultimate state loading condition of the pavement is chosen for all the cases. Figure 10A shows the maximum shear stress profile for the deterministic analysis where the plastic zones are fully developed under the edges of the pavement and the highest magnitude of the maximum shear stress is found to be 36.22 kPa. The plastic zones are symmetric and extended to the bottom of the embankment soil; whereas, the maximum shear stress contours remain no longer symmetric when the spatial variability of the soil properties is considered. Figures 10B,C represent the maximum shear stress profile for the δ x /B δ y /B 0.5, δ z / B 0.5 and δ x /B δ y /B 2, δ z /B 2, respectively with the constant values of COV c 25%, COV tanϕ 20% and L op 3 m.
It is observed that the plastic zones are not fully developed under edges of the pavement for δ x /B δ y /B 0.5, δ z /B 0.5 which may be due to the presence of higher values of shear strength properties of the soil under the pavement; whereas, they are almost developed in case of δ x /B δ y /B 2, δ z /B 2. However, in both the cases they are extended to the mid depth of the embankment. The highest values of maximum shear stress reduced from 45.28 kPa for δ x /B δ y /B 0.5, δ z /B 0.5 to 42.06 kPa for δ x /B δ y /B 2, δ z /B 2 and both the values are comparatively higher than the deterministic value. As per the design point of view, the pavement is considered as unserviceable when the applied stress on the pavement (q app ) is equal or greater than the allowable bearing capacity (q all ) of the pavement. Thus, the limit state of collapse of the pavement can be expressed as q all ≤ q app . In the present situation, q app is the deterministic bearing capacity and q all is the mean probabilistic bearing capacity. Since the assumed distributions for c and tanϕ are lognormal, then the q all is most likely to be lognormally distributed. So, the failure probability of the pavement can be expressed by Eq. 10.
In the above mathematical expression Φ(·) represents the cumulative normal distribution and FOS denotes the factor of safety considered for the design of the pavement. The assumption of considering allowable bearing capacity as lognormally distributed is further assured by performing Kolmogorov-Smirnov (K-S) "goodness of fit" test which is illustrated in Figure 11. The K-S test compares the variation between the actual cumulative frequency and the cumulative distribution function of allowable bearing capacity with assumed theoretical lognormal distribution for L op 3 m, COV tanϕ 20%, COV c 25%, δ x /B δ y /B 2.0 and δ z /B 0.5. The actual frequency of allowable bearing capacity shows a reasonable resemblance with the lognormal fit. Figures 12A,B indicate that the failure probability of the system decreases with the increase in FOS irrespective of horizontal SOF as well as vertical SOF, which is quite obvious. Figure 12A illustrates the variation of failure probability (p f ) for different values of FOS and horizontal SOF with the constant values of L op 3 m, COV tanϕ 20%, COV c 25%, and δ z /B 0.5. For a critical value of FOS 1, the probability of failure reduces with the increase in δ x /B δ y /B. In contrast, it increases with the increasing values of δ x /B δ y /B for higher factor of safety. This can be attributed as while the FOS 1, for small values of δ x /B δ y /B, the COV q is comparatively smaller than the larger values of δ x /B δ y /B. Therefore, the µ q tends to drop below the limiting value and the failure probability increases. However, with the increase in δ x /B δ y /B, the failure probability decreases due to the increase in stability for FOS 1. For higher values of FOS, the µ q is found to be pushed above the limiting value which reduces the failure probability for small values of δ x /B δ y /B. However, the increase in δ x /B δ y /B causes increasing instability which in turn increases the failure probability of the system for higher FOS.
Similar type of response has been observed in Figure 12B which represents the variation of failure probability (p f ) for different values of FOS and vertical SOF with the constant values of L op 3 m, COV tanϕ 20%, COV c 25%, and δ x /B δ y /B 2. However, in this figure, the failure probability is found to be decreasing with the increasing δ z /B for FOS 1 and 1.1. For FOS greater than 1.1, the failure probability is increasing with the increase in δ z /B. Beyond a certain scale of fluctuation, the failure probability is observed to be almost insensitive with the scale of fluctuation for both the figures. The similar kind of trend has been observed by Griffiths and Fenton (2004), Chenari and Alaie (2015) and, Halder and Chakraborty (2020). Figures 12A,B also compare the failure probability between the consideration of spatial variability and without consideration of spatial variability for different values of FOS. It is evident from the figures that the failure probability without consideration of spatial variability does not depend upon the scale of fluctuation. It is also observed in Figure 12A that for FOS 1, the failure probability without considering the spatial variability is lower than that of considering the spatial variability. For FOS 1.1, the probability of failure for smaller values of δ x /B δ y /B is found to be lower as compared to that of without consideration of spatial variability; whereas, for higher values of δ x /B δ y /B, it becomes higher than that of without considering the spatial variability. However, for FOS greater than 1.1, the failure probability without considering the spatial variability is found to be quite higher as compared to that of considering the spatial variability. For an instance, p f drops from 86.99% with the consideration of spatial variability (δ x /B δ y / B 2 and δ z /B 0.5) to 60.4% without considering the spatial variability for FOS 1. For FOS 1.1, the failure probability for δ x / B δ y /B 1.5 and δ z /B 0.5 becomes 37.4% which is lower as compared to 43.9% (p f corresponding to without considering the spatial variability); whereas, it increases to 50.2% for δ x /B δ y /B 10 and δ z /B 0.5. However, for FOS 1.2, p f increases from 14.4% with the consideration of spatial variability (δ x /B δ y /B 4 and δ z /B 0.5) to 29.7% without undertaking the spatial variability. Similar kind of responses are observed in Figure 12B. However, unlike δ x /B δ y /B, in case of δ z /B, the probabilities of failure without considering the spatial variability are found to be higher for FOS greater than 1, as compared to that of considering the spatial variability. As an example, p f decreases from 73.8% with the consideration of spatial variability (δ x /B δ y /B 2 and δ z /B 2) to 60.4% without considering the spatial variability for FOS 1; whereas, for FOS 1.2, it increases from 10.9% with the consideration of spatial variability (δ x /B δ y /B 2 and δ z / B 8) to 29.7% without undertaking the spatial variability.

CONCLUSIONS
The present paper presents the three-dimensional probabilistic bearing capacity-settlement response of the pavement located on the crest of the embankment where the embankment soil is considered to be fibre-reinforced. The deterministic as well as the probabilistic studies are executed. The probabilistic parameters considered for the random field generation are taken from the previous studies. In probabilistic analyses, the effects of soil spatial variability, randomness of soil properties and different out of the plane lengths have been studied. The conclusions drawn from the study are as follows: 1. The mean bearing capacity of the pavement is found to be decreased first with increasing values of the out of the plane length of the embankment (L op 0.5-1 m), and after that, it starts to increase with the increase in L op (L op 1-3 m).
Beyond L op 3 m, there is no significant change in it. In contrast to that, the COV q is observed to be decreasing with the increasing value of L op . It signifies that the consideration of different out of the plane lengths with spatial variability is quite important as it makes the problem more realistic. 2. The mean bearing capacity is found to be reducing with the increment in randomness in the soil friction angle (COV tanϕ ); whereas, the COV q increases with the increase in randomness in the soil friction angle which is quite obvious as it leads to the lower strength paths for failure. 3. The mean and COV of bearing capacity of the pavement are turned up to be decreasing and increasing, respectively with the increasing values of horizontal SOF. Unlike the horizontal SOFs, the mean and COV of bearing capacity both are found to be increasing with the increasing vertical SOFs. 4. The plastic zones under the pavement edges are found to be fully developed for homogeneous soil; whereas, they are partially developed for spatially variable soil. The highest values of maximum shear stress for spatially varied soil are found to be higher than that of soil with uniform strength FIGURE 11 | Comparison between actual distribution and the assumed theoretical lognormal distribution of allowable bearing capacity for L op 3 m, COV tanϕ 20%, COV c 25%, δ x /B δ y /B 2 and δ z /B 0.5. properties which implicates the importance of using spatial variation of the soil properties. 5. The failure probability of the pavement is found to be decreasing with the increase in scale of fluctuation for lower values of FOS; whereas, it is observed to be increasing with the increasing values of scale of fluctuation for higher values of FOS. 6. The failure probability of the pavement is observed to be 29.7% for FOS 1.2 when the spatial variation of soil properties is not considered; whereas, it is found to be reduced to 9.73% for FOS 1.2, δ x /B δ y /B 2, and δ z /B 2. This further signifies the importance of considering spatial variability of the soil shear strength properties.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
KC generated the MATLAB and FISH codes to carry out the probabilistic study and prepared the initial draft. DC proposed the project, supervised all the works and revised the initial draft to give the final shape of the paper.