Effect of soil anisotropy and variability on the stability of undrained soil slope

Soil is a naturally heterogeneous material and can show significant spatial variation in strength and other properties. For silty and clayey soils, these variations are often more pronounced. Despite such variation, many past studies considered these soils as homogeneous and only considered a single set of soil parameters. This may lead to underestimation of the failure potential of geo-structure such as natural slopes, water retaining dams, retaining walls, etc. A finite element method considering soil variability should be an ideal tool to investigate the behaviour of these soils. This study adopted a 2D random finite element method to evaluate the effect of such variability on slope stability. The spatial variability was implemented by using the coefficient of variation (COV) and the spatial correlation length (θ) for cohesion. It was found that the soil slope with higher COV would have a higher chance of failure, whereas the soil slope with less COV might not show any failure. In addition, the soil with a higher θ, in general, show less potential of failure. In the literature, most studies considered an isotropic condition for the soil, i.e., θ in x and y directions are the same θ x = θ y , which is not realistic. Therefore, the soil anisotropy (i.e., θ x ≠ θ y ) was considered carefully in this study. It was found that the probability of failure for anisotropic soil might be significantly higher than the isotropic soil.


Introduction
Soil slopes, i.e., cuttings or embankments, are an integral part of our transport network. Since the variability of material quality is common in soils, especially for silty or clayey soils (Karim et al., 2011;Devkota et al., 2022;Karim et al., 2022), the uncertainty on the stability of these slopes is an inherent challenge (Griffiths, 1982;Huang et al., 2006;Gould et al., 2011;Kasama and Whittle, 2011;Karim et al., 2021). Further, the mechanics for slope stability analysis that involves the conventional limit state approach uses simplifications and approximations, which adds to the uncertainty of their design. Thus designers often use a highly conservative estimation of strength parameters (Griffiths and Lane, 1999). The assumption of homogenous soil is also common in geotechnical designs, even though the strength parameters may vary spatially (Idriss et al., 1978;Huang et al., 2013;Rahman and Nguyen, 2013). Therefore, an alternative approach that allows natural soil characteristics of spatial variability is needed to quantify the related risk of failure.
In recent years, an advanced numerical technique that uses a combination of the finite element method (FEM) and random field model (RFM) has been gaining popularity and is capable of capturing the soil characteristics of spatial variability in soil properties (Griffiths and Lane, 1999;Griffiths and Fenton, 2001;Fenton and Griffiths, 2003;Fenton and Griffiths, 2008;Rahman and Nguyen, 2012;Huang et al., 2013;Li et al., 2019). Griffiths and Fenton, (2001) adopted a technique called local average subdivision (LAS) for generating a random field for undrained shear strength. This technique was first introduced by Fenton (1990), and then it was revisited by others (Fenton and Vanmarcke, 1990;Fenton and Griffiths, 2007) based on the probability density and spatial correlation functions for soil parameters. The main advantage of this approach is the capability of capturing the complex behaviour of soil with high spatial variability.
In this study, the FEM with a random field generator is used for slope stability analysis. This method is also referred to as the random finite element method (RFEM), as proposed by Griffiths and Fenton, (2001). RFEM has been applied in many different areas of geotechnical engineering, such as bearing capacity under foundation, slope stability, pile foundation, etc. (Smith and Griffiths, 1998;Griffiths and Lane, 1999;Fenton and Griffiths, 2008;Jiang et al., 2022a;Jiang et al., 2022b;Nguyen et al., 2022;Shu et al., 2023). It allows the soil properties to be changed spatially but still correlated with the neighbouring soils (Fenton and Vanmarcke, 1990;Fenton and Griffiths, 2007). Griffiths and Fenton, (2001) proposed correlation length (θ) as a measure for investigating the correlations between neighbouring soils. Based on the correlation length and soil variability, the random field for soil parameters can be generated. One of the most efficient random field generation techniques is the local average subdivision (LAS), which was introduced in Vanmarcke (1984). This technique considers a global average for a parent cell. Then, the cell is subdivided multiple times. The local average in each child cell is recalculated. This technique ensures a reasonable transition in soil parameters between the neighbouring cells. It should be noted that there have been different studies considering other techniques for generating a random field. One notable approach, the numerical limit analysis (Kasama and Whittle, 2011) (NLA) is based on the numerical formulations of upper and lower bound limit analyses for rigid perfectly plastic materials, using finite element discretisation and linear (Sloan, 1988;Sloan and Kleeman, 1995) or non-linear programming methods (Lyamin and Sloan, 2002a;Lyamin and Sloan, 2002b). According to Kasama and Whittle (2011), both upper and lower bound analysis treats the soil elements as three-node triangular elements. By applying the upper and lower bound limits, the calculated failure load was covered into a range. However, the NLA technique may produce results that do not align well with the theoretical data, whereas the LAS technique can produce the best-fit results with the theoretical values. Therefore, this research considered the LAS technique for generating a random field.
Most RFEM studies in slope stability considered an isotropic condition in which the correlations between neighbouring soils in lateral and vertical directions are the same (Griffiths and Fenton, 2004;Griffiths et al., 2015;Kasama and Whittle, 2016;Zhu et al., 2017). However, an anisotropic condition is more common in nature. Soils in their natural states deposit in layers, and their strength and other properties can be different in vertical and lateral directions, i.e., anisotropic conditions prevail. Most previous studies often found the correlation lengths in the horizontal direction were higher than in the vertical direction (DeGroot and Baecher, 1993;Vessia et al., 2009;Rahman and Nguyen, 2013).
There have been some studies that incorporated soil anisotropy (Rahman and Nguyen, 2013;Nguyen and Rahman, 2015;Pieczyńska-Kozłowska et al., 2015;Li et al., 2016;Zaskórski et al., 2017) in bearing capacity problems. Some studies were done for slope stability. Akbas and Huvaj (2015) only evaluated different combinations of correlation length with a single value of the coefficient of variation (COV) of 0.2, which represented slightly varied soil. Liu et al. (2018) investigated the stability of 3D slope with one combination of correlation length (θ x = 2H and θ y = 0.4H, where H is the height of the slope) and COV of 0.3 for undrained cohesion. The studies using COV less than 1.0 may not be able to capture the randomness of highly varied materials such as tailings, etc. Other studies (Huang et al., 2013;Li et al., 2019) examined the finite element strength reduction analysis in RFEM; however, the ranges of COV and correlation are not wide enough to explore different possibilities of a factor of safety. Therefore, the combined effect of soil anisotropy and variability on the probability of slope failure is not fully understood, mostly due to a limited number of past studies incorporating anisotropy and limited data sets available. Therefore, this research adopts the RFEM approach from Griffiths' group (Griffiths and Lane, 1999;Griffiths and Fenton, 2001) with LAS to consider the soil variability, simulates both isotropic and anisotropic conditions by combining a wide range of coefficients of variation (COV) of soil properties and correlation lengths in the vertical and horizontal directions and then investigate the combined effects of correlation length and soil anisotropy on slope stability. This study produces a comprehensive data set of anisotropic conditions for correlation lengths for the main investigation.

A probabilistic method for slope stability
The stability of a soil slope in the undrained condition is often analysed using a set of soil properties/parameters consisting of the undrained shear strength (c u ), the undrained friction angle (ϕ u ), and the saturated unit weight (γ sat ). In the deterministic approach, these parameters are assumed to be constant, i.e., a homogenous soil. The friction angle, ϕ u , for silty or clayey soils can be taken as zero. For An illustration of the slope geometry in this study.

Frontiers in Built Environment
frontiersin.org 02 simplicity, γ sat is also held constant at 20 kN/m 3 . Other important parameters for the slope are the slope inclination (ß), the height of slope (H) and the depth of foundation indicated by n d . An illustration of the slope used for analysis is shown in Figure 1. The soil strength parameters are reduced during the analysis to capture the failure of the slope.
The undrained cohesion (c u ) is allowed to vary spatially, and its effect on slope stability is examined. To compare with the deterministic approach, c u is generalised and expressed in terms of a dimensionless parameter (C) by normalising against the multiplication of unit weight and slope height.
Two different approaches are adopted to generate the spatial variation of C in the slope geometry, namely, the single random variable approach and the random variable approach considering spatial correlation length. In the single random variable approach, a lognormal distribution function for C was adopted to avoid the unreasonable negative value of C and to ensure positive shear strength by the random field generator. After Griffiths and Fenton (2004), the probability density function of C is expressed as a log-normal distribution function (see Figure 2): where μ lnC and σ lnC are the mean and standard deviation of the natural logarithm of C. μ lnC and σ lnC are determined by the dimensionless coefficient of variation (COV), which is defined as the ratio μ C and σ C as shown below.
According to Eq. 1, the COV of the dimensionless parameter (C) was equal to COV cu (COV cu σ cu /μ cu ), as γ sat and H were kept constant in this study. The distribution of ln(C) is then characterised by: To continue with the calculation of the probabilistic failure, the chance of failure has to be defined by the probability of random C less than the characteristic value of C. For this, a random value of C from the probability function of Eq. 2 was assigned for the entire slope. Therefore, the slope is still homogeneous. According to Griffiths and Fenton (2004), the relationship between C and the factor of safety (FOS) of the homogeneous slope can be presented in Table 1. The data in Table 1 shows a linear relationship (plotted in Figure 3) between C and FOS with a very high coefficient of determination (R 2~1 ).
Based on the data in Table 1 and Figure 3, the characteristic value of C for a slope with FOS of 1.00 is 0.17. Hence, the chance of failure of the slope in this probabilistic study can be expressed as: In the second analysis with the spatial correlation, one of the important parameters will be the spatial correlation length (θ lnC ) (Griffiths and Fenton, 2001). When C is log-normally distributed, its logarithm yields an underlying Gaussian field. θ lnC , which is the measure of this Gaussian field, ensures that the neighbouring soils in the slope are still correlated. A normalised dimensionless measure of θ lnC (Θ lnC ) can be presented by the ratio of θ lnC and H. The representative value of θ lnC for clay has not been clearly defined in the literature (Griffiths and Fenton, 2001), especially in the horizontal directions. This is because the soil normally deposits into layers and the higher variation in the vertical direction is more likely to be observed. This analysis is presented in the following section.
It should be noted that the Monte Carlo simulation was adopted in this study. This method uses the estimation of randomly stochastic soil property based on its distribution type, mean and standard deviation. The Monte Carlo simulation also uses the correlation function or spectral density function, which is also a characteristic of spatial variability for cohesion. Each test considered 1,000 Monte Carlo simulations.

Effect of soil variability and anisotropy on slope stability
As mentioned previously, the elastoplastic FEM coupling with the probabilistic approach is an effective tool that can better capture the effect of variability of soil in the slope stability problem (Akbas and Huvaj;Griffiths and Fenton, 2004;Vessia et al., 2009;Huang et al., 2013;Rahman and Nguyen, 2013;Pieczyńska-Kozłowska et al., 2015). The discretised geometry used in this study is presented in Figure 4 and the LAS technique is adopted for generating the random field due to its efficiency and accuracy. The mean undrained cohesion is chosen as 100kPa, and the standard deviations for the undrained cohesion are 25kPa, 50kPa, 100kPa, 200kPa, 400kPa and 800 kPa respectively.
Each mesh in Figure 4 has a corresponding value of C throughout the generating process in the LAS technique. The Log-normal distribution function for C.

Frontiers in Built Environment
frontiersin.org correlation function between two distant points is defined by the Markovian function (Fenton and Griffiths, 2008) as where X ij and Y ij are the distances between two random points in horizontal and vertical directions respectively. The correlation field of C will be different with different combinations of θ x andθ y . Figure 5 demonstrates the capacity of different combinations of θ x andθ y to generate the correlation field or correlation coefficient (ρ). Figures 5A-C show the symmetrical distributions of ρ for the isotropic conditions.
• For a lower correlation length (θ x = θ y = 1.0), the highest correlation is achieved when the random points are closer to each other i.e., X ij and Y ij are nearly zero. When the two random points are far away from each other, ρ reduces significantly. • For a higher correlation length (θ x = θ y = 10 or 100), ρ only reduces slightly, when the two random points are far away from each other.
In Figures 5D-I, the distributions of ρ for the anisotropic conditions are shown. When θ x >θ y , ? y reduces faster than ? x in the case of two far random points, and vice versa. It should be noted that when θ x orθ y is very high and approaching infinity, the correlation coefficient, ρ is approaching 1. This means 100% correlation between two random points, which is an optimistic assumption in the design.
The correlation coefficient (ρ) between separated points in anisotropy was much different with different combinations of the correlation lengths. ρ tends to increase to 1 in the x-direction while θ x is continuously increasing up to infinity. The random set of variables obtained from the correlation function is later used to transform the random field of c within the random field with available mean and standard deviation. The lognormal transformation is defined as  (Griffiths and Fenton, 2004

FIGURE 3
The relationship between C and FOS of a homogeneous slope.

FIGURE 4
The discretised geometry of the slope for RFEM study.
Frontiers in Built Environment frontiersin.org where G is a standard lognormal variable, which is defined throughout the LAS process. Then the function used to subdivide the mesh in the LAS process is adopted: where U is indicated as the vector of independent standard normal random variables with mean zero and unit standard deviation. The covariance describing the relationship between the cells can be calculated by the matrix multiplication with a transpose

R E QQ T (9)
S E QG T (10) where R is the covariance between the parent cells, S is the covariance between the parent and child cells and B is the covariance between the child cells. Then the matrices A and L in Eq. 8 can be determined by The study considers a wide range of θ x ,θ y and COV to investigate the combined effect of correlation length and COV on the probability of failure. One combination of θ x ,θ y and COV gives one p f . In total, there are 216 simulations in this study. The details can be found in Table 2.
The slope with a small variation of soil parameters (i.e., low COV of 0.25-0.50 and high θ of 8 m), a smooth random field can be observed (see Figure 6). Figure 6A shows the variation of c u and the deformations of the cells; whereas Figure 6B shows only displacement vectors. It should be noted the darker colour in Figure 6A means a higher value for cohesion. There is no failure surface appearing in this case, as the displacement vectors do not reach the right side of the slope in Figure 6B. As mentioned before, the elements were assigned different C parameters. The failure mechanism will pass through these elements by find the weakest zone based on the assigned parameters. The slope will fail if the deformation or displacement is significant. It noted that this research also considered the shear strength reduction technique for the slope stability analysis (Griffiths and Lane, 1999).
On the other hand, if the soils have high variation (i.e., high COV) and small correlation (i.e., small θ), the cohesion field is ragged. In this case, cohesion values of neighbouring soils are mostly random at the final state of the slope analysis as shown in Figure 7A, and the probability of failure is 100%. The failure line appears on the slope and follows the weakest path traversing the weakest points in the slope geometry. The directions of the displacement vectors also follow the same weakest path as shown in Figure 7.
The traditional approach for slope stability analysis often assumes a deterministic value for soil strength as mentioned previously. This means that the soil parameters have no variation, which is similar to the case presented in Figure 6. In such a case, the chance of failure is very low. However, the assumption of slight variation for a parameter is not always true for most soils, especially mining materials. A big failure surface appears when the variation increases (see Figure 7). Hence, it may lead to underestimating the chance of slope failure in traditional slope analysis.   Although acknowledging the variation of soil properties is important for engineering design, the current practice is still based on a deterministic approach. The probability of slope failure (p f ) is found to depend on COV and θ as shown in Figure 8. With a small COV, i.e. 0.50, the p f is relatively small. The p f value is approximately zero at small θ i.e. 0.20 m and increases with an increase in θ. On the other hand, when COV is higher (from 1.00 to 4.00), the trend of p f seems reversed. At COV = 4.0, the p f value at a small θ of 0.20 m is approximately 100%. This value decreases when θ increases. This demonstrates that the soil with higher soil variability (high COV and low θ) will have a much higher chance of slope failure. Figure 8 also shows the difference between isotropic and anisotropic θ. In some cases, considering an anisotropic condition, which is more realistic, may give a higher value of p f . Using isotropic θ, which is mostly considered in the literature, may not be conservative in the case of highly varied soil and the potential failure may not be accurately determined. Figure 8E shows the surface map of the 3D distribution of p f based on different combinations of COV and θ. It can be clearly seen that the p f surface bends upward when θ increases for low COV, i.e. 0.5. On the other hand, the p f surface bends upward when θ increases for COV > 1.0. This clearly demonstrates that low COV and high θ may result in higher p f , whereas high COV and low θ may result in higher p f . According to Eq. 6, the correlation between the points is calculated by θ. The highest correlation with small θ can be achieved, when the points are closer to each other. In that case, the FOS values from Monte Carlo simulations will be mostly higher than the FOS values and p f is closer to 0. When the two random points are far away, the correlation can be reduced significantly. A similar observation was found in the study of the isotropic slope by Griffiths and Fenton (2004).
Furthermore, p f is further investigated by considering COV and dimensionless θ (Θ) in Figure 9. The effect of COV ranges from 0.5 to 8.0 on the failure potential of the slope in isotropic conditions. The probability of slope failure, p f can increase from less than 20% (COV = 0.5) to near 100% (COV = 8.0). This indicates that the soil with high variability can have a much greater p f compared with the soil with less variability. It should be noted that the x-axis in Figure 9 presents the dimensionless correlation length, Θ, which is the ratio

FIGURE 9
Effect of COV and θ on the probability of failure in isotropic conditions.

FIGURE 10
Effect of COV and θ on the probability of failure in anisotropic conditions: (A) θ y 0.5 m, (B) θ y 1.0 m and (C) θ y 4.0 m.

Frontiers in Built Environment
frontiersin.org of spatial correlation length to the slope height. It can be observed that p f increases with Θ, when COV is less than 1. But this trend changes when COV is more than 1. So, p f of a soil slope is not affected by a single factor, but by the combination of soil variability and correlation. This should be carefully considered in slope design.
As mentioned previously, this study also considers the effect of soil anisotropy on the probability of slope failure. Figure 10 presents the combined effect of COV and anisotropic θ. In Figure 10A (θ y = 0.5), the p f values are almost 1.0, i.e. 100% chance of failure for COV ≥ 2.0. However, when θ y increases, the p f values start diverging at high COV (see Figures  10B, C). This observation aligns with the previous discussions that the more correlation length is, the less chance of failure is.
In addition, Figure 10 presents the combined effect of the degree of anisotropy (i.e., ratio of θ x and θ y ) and COV on p f . Interestingly, the lines representing different degrees of anisotropy crossover each other in Figure 9. Before the crossover, the lower θ x /θ y ratio shows a lower p f . After the crossover, the lower θ x /θ y ratio shows a higher p f . The location of the crossover is highly dependent on the combination of correlation length and COV. For instance, the crossover for θ y = 0.5 occurs at COV of 0.5, whereas the crossovers for θ y = 1.0 and 4.0 occur at COV of 0.6 and 0.75 respectively. So, such a reverse effect is clearly dependent on the θ y and COV. This is an interesting finding that may help in the future study of risk assessment of slope stability. However, it should be aware that other soils with more varying parameter sets (e.g., varying friction angle, stiffness, unit weight, etc.) may be different.

Conclusion
The study introduces a more reliable approach for the quantification of risk in slope design. The combination of a random field generator and FEM (RFEM) considers the variability of soil parameters, which is more realistic than the common deterministic approach. The probability of slope failure is found to be greatly dependent on the statistical variation parameters, namely, COV and θ. This study produces a wide range of different scenarios that can happen in the real design and outputs a set of different p f values for a slope while the traditional approach assumes a p f value. Hence, this approach can help to minimise the chance of overestimation or underestimation of slope failure. The effect of anisotropy, which is often neglected in the current engineering practice, is also proposed in this research. In some cases, the p f values of anisotropic θ are much higher than that of isotropic θ, which helps to prevent possible design failure.
Interestingly, this was found that the combination of a low COV and a high θ may result in a higher p f , and vice versa. The knowledge from this RFEM study can be adopted in any future research in the risk assessment of geo-structures such as designing a slope of soil with high variability such as tailings. In those cases, the practitioners may consider the additional factor of safety or some possible soil improvement technique to avoid potential failure.
It should be noted that for most natural soils, a reasonable COV range is from 0.1 to 0.5. However, for the materials having very high variability such as tailings, the COV value should be higher.

Data availability statement
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.