Abstract
This paper concentrates on the fact that natural materials are variable, and that representation of this variability appears crucial to getting a realistic understanding of certain geotechnical problems. In this case, finite elements have been used to assess the influence of heterogeneity on the stability of a clay slope, characterized by a spatially varying undrained shear strength cu. For low to intermediate values of the coefficient of variation (e.g. 0.1–0.3), as are often observed in practice, pointwize variability can be approximated by a normal distribution. Data on spatial variability are more scarce, although vertical scales of fluctuation are usually much smaller than horizontal scales of fluctuation, and will often be small relative to the height of a slope. The results show that, for a given factor of safety, reliability is greatest for cu constant with depth. For cu increasing linearly with depth, reliability decreases due to the greater range of possible rupture surfaces: hence, in conventional (deterministic) design, higher factors of safety (or lower characteristic values) would be needed. In this study, horizontal anisotropy (of the heterogeneity) causes modest changes in reliability. However, problems are identified in which anisotropy is likely to have a significant influence.
Soil structures have traditionally been analysed by subdividing into zones (or layers) of uniform material. However, it is well known that, even in so-called ‘uniform’ deposits, spatial variations in material properties do exist (i.e. heterogeneity). Hence, at any point in a uniform layer of clay, for example, the undrained shear strength will be some function of a statistical distribution, rather than some fixed value which is constant throughout the layer.
The conventional (deterministic) approach is to adopt characteristic (e.g. mean) property values for each layer. This leads to a single analysis and, for stability assessments, leads to a single factor of safety. In contrast, the stochastic approach takes account of all properties within each layer. This requires multiple analyses, as part of a Monte Carlo simulation process, and leads to an alternative, more meaningful, definition of stability: i.e. reliability, the probability that failure will not occur.
This paper studies the influence of heterogeneity on slope stability, by linking random field theory, for generating property distributions, with finite elements, for assessing stability. This approach was used by Hicks & Onisiphorou (2000) and Onisiphorou (2000) for assessing the liquefaction potential of an underwater sandfill slope. However, a more familiar problem is considered here: specifically, an undrained clay slope, characterized by the undrained shear strength cu. This problem has been studied previously, by Paice & Griffiths (1997) and Griffiths & Fenton (2000), who based their parametric studies on an isotropic spatial variation of cu relative to a mean which was constant with depth.
In this investigation, the influence of mean undrained shear strength varying with depth has been considered, as has anisotropic heterogeneity, which is apparent in all geomaterials due to the process of deposition. An important aspect of this study has been the use of statistical parameters based on laboratory and field observation.
Summary of stochastic process
Figure 1 illustrates the statistical properties of an undrained clay, characterized by the undrained shear strength cu. Firstly, Figure 1(a) shows cu varying with depth in a uniform deposit. Conventional design would make use of either the mean undrained shear strength, which, for this example, is constant with depth, or some other characteristic value. In contrast, stochastic design makes use of all data, as illustrated in Figure 1(b), in which the data are expressed in the form of a probability density function (pdf), characterized by the mean undrained shear strength, μ=μcu, and by the standard deviation from that mean, σ=σcu. These point statistics lead to the coefficient of variation, COV=COVcu, in which
Statistical measures of cu variability. (a) cu varying with depth, (b) probability density function of cu.
The pdf in Figure 1(b) follows a normal distribution, although other distributions are sometimes used: e.g. uniform, lognormal, exponential and beta (see Lee et al. 1983). Furthermore, the definitions of mean and standard deviation can be extended so that they arefunctions of depth z. For example, Lumb (1966) considered a linear trend, i.e. for an undrained clay,
in which
and where u is a standardized random variable of zero mean and unit standard deviation, which varies according to the type of pdf: for a normal distribution, this is given by,
In equation (3), for mean cu to be constant with depth (as in Figure 1(a)), βμ=0, while αμ=ασ=0 causes both μ and σ to increase linearly from values of zero at z=0: in this case, COV is constant with depth and given by βσ/βμ.
Equations (1)–(4) define pointwize variability. Extension to spatial variability requires a further statistical parameter, θ, defining the degree of spatial correlation, or distance beyond which there is minimal correlation, e.g. Vanmarcke (1983). Figure 1(a) shows that θ is a function of the distance between adjacent weak or strong zones. Hence, for a weakly correlated material θ is small, while, for a strongly correlated material, θ is large.
Once the statistics (μ, σ, θ) for a material are known, a numerical simulation of the material property distribution (i.e. a random field) can be generated. In this paper, this has been done using local average subdivision (LAS) (Fenton & Vanmarcke 1990), which generates correlated local averages, Z, based on a Gaussian pdf (zero mean and unit variance) and spatial correlation function ρ. For an isotropic random field, i.e. θ equal in all directions, and for a square domain, the Gauss–Markov correlation function is given by,
in which τ is the lag vector. The random field is then transformed from a Gaussian field into that defined by the actual pdf, by using μ(z) and σ(z) from equation (3). For a normal distribution, the transformation is,
in which, for 2-D problems, x=(x z)T, and Z(x) is the local average for a random field cell with centroid x. LAS generates a square random field by uniformly subdividing a square domain into smaller square cells, each cell having a unique local average that is correlated with surrounding cells.
For a given set of statistics, there are an infinite number of possible random fields: although statistically similar, each will yield a different solution, depending on the spatial distribution of material properties. Hence, stochastic analysis involves the repeated analysis of realizations as part of a Monte Carlo simulation process. The results may then be expressed as a pdf, enabling stability to be expressed, either in terms of risk (the probability of failure), or reliability (the probability of failure not occurring). In any such simulation, enough realizations are therefore needed to ensure convergence of the output statistics.
Evidence for statistical parameters
The present investigation has similarities with previous studies by Paice & Griffiths (1997) and Griffiths & Fenton (2000), who examined the influence of isotropic spatial variability on factors of safety for undrained clay slopes (of height H), for μ constant with depth (i.e. the Taylor (1937) problem). Based on a lognormal property distribution for cu, the following parameter ranges were considered: 0.2≤COV≤1.0 and H/40≤θ≤∞ (Paice & Griffiths 1997), and 0.125≤COV≤8.0 and H/2≤θ≤∞ (Griffiths & Fenton 2000).
This investigation focuses on statistical evidence from laboratory and field test data. For example, Lumb (1966) used the χ2 test, and data from vane tests and unconfined compression tests, to demonstrate that cu variability follows a normal distribution: similar findings for cu have been found by others, e.g. Chiasson et al. (1995) and Hooper & Butler (1966). Lumb (1966) further demonstrated that normal distributions were suitable for other soil properties (e.g. effective cohesion and friction angle), while Onisiphorou (2000) analysed data from 71 CPTs on sands to demonstrate the normality of state parameter variability (Been & Jefferies 1985). The applicability of the normal distribution, for soil properties in general, was also supported by Lee et al. (1983), and has usually been used by investigators when defining cu variability.
For stochastic analysis, the normal distribution is conceptually simpler than most other distributions and problems associated with skewness do not arise. For example, for higher COVs, the increased skewness of lognormal distributions can lead to an increased difficulty in interpreting results, as well as a substantial increase in the number of realizations needed for convergence of the output statistics. A disadvantage of the normal distribution is the possibility of negative property values; a problem not encountered with the lognormal distribution, though easily rectified for low and intermediate values of COV, i.e. by stipulating a lower bound for cu. For a normal pdf, 99.73% of data lie in the range, μ±3σ: hence, for COV≤0.33, less than 0.135% of data will need to be truncated.
Lee et al. (1983) and Phoon & Kulhawy (1999) suggested a likely range for COV (for cu) of 0.1–0.5. However, other published values, based on in-situ testing, include: 0.18 and 0.16, for marine and London clays, respectively (Lumb 1966); 0.21 and 0.23, for two marine clays (Soulié et al. 1990); 0.18 to 0.30, for New York City clay (Asaoka & Grivas 1982); 0.32 for varved clay (Vanmarcke 1977b); while a value around 0.2 may be inferred from the data of Chiasson et al. (1995). Also, Hooper & Butler (1966) reported values ranging from 0.11 to 0.33, based on an extensive laboratory investigation of London clay. In this case, the results could be divided into two distinct groups, according only to the method of sampling. Both groups showed similar trends in mean cu, but average COVs of 0.14 and 0.28, therefore raising the possibility that higher values of COV may sometimes be an indication of sampling technique, rather than actual (inherent) material variability. Furthermore, if a deposit shows a mean cu increasing with depth, as is often the case, then this trend should be removed for computing COV—otherwise, COV may be overestimated and, if σ increases with depth, a skewness imposed on the pdf.
Spatial correlation data are more scarce, and interpretations more subjective, due partly to the adopted definition of correlation (e.g. scale of fluctuation, autocorrelation distance), and partly to the evaluation method. For example, Popescu et al. (1995) computed vertical scales of fluctuation of 1.8 m and 0.84 m from the same NSPT data for a sandy soil at Akita Harbour, based on the correlation function and variance function methods, respectively. Other methods have been proposed by Wickremesinghe & Campanella (1993), who reported a scale of fluctuation of 0.3 m for cone bearing in sand, and Soulié et al. (1990) and Chiasson et al. (1995), who computed autocorrelation distances for cu of 3 m and 2 m, for marine and sensitive clays, respectively. Onisiphorou (2000) used Wickremesinghe & Campanella’s (1993) method in computing scales of fluctuation, for state parameter, of around 0.7 m for a 26 m-high hydraulic sandfill slope, while Vanmarcke (1977b) reported a θ for cu of 5 m for a deep clay deposit, based on Vanmarcke (1977a). A further uncertainty may lie in the nature of spatial variations being measured: i.e. larger correlation values may sometimes refer to spatial variations linked to layering, rather than to variations within layers themselves, as has been considered in this paper. Furthermore, non-removal of mean trends can also lead to overestimation.
Data for correlations in the horizontal plane are rare. However, available evidence strongly points to a much higher degree of correlation than in the vertical direction, due to the process of deposition. For example, Popescu (1995) and Popescu et al. (1997) compared closely-spaced CPT profiles taken from a hydraulic sandfill and found scales of fluctuation in the horizontal direction to be, on average, 12 times greater than in the vertical direction. Similar findings have been reported for the degree of anisotropy of clays: e.g. 9 (Vanmarcke 1977b), 10 (Soulié et al. 1990), and 13 (Phoon & Kulhawy 1999).
Numerical simulation
Figure 2 summarizes the problem geometry and boundary conditions for a 45° clay slope, of height H=5 m, founded on a firm base. Its stability has been assessed using finite elements and the method outlined by Smith & Griffiths (1998): hence, the in-situ stresses in the slope have been generated by applying gravitational forces to the finite element mesh, based on the soil's unit weight, γ=20 kN/m3 (cf. the loading of a model slope in a centrifuge). The clay has been modelled using an elastic, perfectly plastic, Tresca soil model.
Problem geometry, boundary conditions and mesh discretization.
Figure 2 shows that the problem domain has been discretized using 610, 8-node, quadrilateral elements, each with 2×2 Gaussian integration, and that, apart from along the sloping boundary where elements have been distorted to form triangles, the mesh comprises square elements of size 0.25 m (H/20). This high level of discretization is necessary for the accurate modelling of random fields with small scales of fluctuation. However, optimal use of elements has been achieved by assigning random field (cell) values to element sampling points, rather than to elements themselves: i.e. a group of four random field cells occupies the same area as one finite element, and each realization therefore uses 2440 values of cu (i.e. 4×610). This increase in efficiency is possible due to the accuracy of the 2×2 Gaussian sampling point locations for 8-node elements (Onisiphorou & Hicks 2001).
All property distributions have been obtained from random fields of cell size 0.125 m (H/40), which have been generated for a square domain (larger in area than the problem domain). The sequence of assigning properties for an isotropic distribution is as follows:
-
for the given value of θ, and using equation (5), generate an isotropic Gaussian random field Z(x), i.e. mean zero and variance unity;
-
transfer the cell values directly onto the sampling points in the problem domain;
-
transform the sampling point values into cu, using μ(z), σ(z), and equations (6) and (3).
Anisotropic distributions are obtained in a similar manner:
-
generate an isotropic Gaussian random field based on θ=θh, the scale of fluctuation in the horizontal direction;
-
distort the random field in the vertical direction by compressing the cells by a factor ξ=θh/θv, the degree of anisotropy;
-
average cell values in the vertical direction so that new square cells are produced; hence, the modified random field will be rectangular—in the vertical direction it will be ξ times shorter than in the horizontal direction and contain ξ times fewer square cells (note that the distortion process has been made simpler by considering only cases in which ξ is a whole number);
-
transfer cell values (as above);
-
transform sampling point values (as above).
Figures 3 and 4 show typical random fields of cu, for μ constant with depth and increasing linearly with depth, respectively. In these figures, light colours represent lowvalues of cu, and dark colours, high values of cu. Using laboratory and field evidence as a guide, a relatively narrow range of statistics has been considered: specifically, for scales of fluctuation in the vertical direction, 0.25≤θv≤2.5 metres, or H/20≤θv≤H/2, while, for θv=1.0 m (=H/5), the range of anisotropy considered is, 1.0≤ξ≤∞ (note that cu distributions for ξ=∞ have been generated using a 1-D LAS algorithm). For all stochastic analyses, COV=0.3: although this is towards the upper end of probable values, it is also small enough to ensure that truncation of the normal distribution has minimal effect (note that a lower limit of 1 kPa has been set for cu). The elastic soil properties are: Young's modulus, 105 kPa, and Poisson's ratio, 0.3.
Typical random fields of cu for μ constant with depth. (a) θv=1.0 m=H/5; ξ=1, (b) θv=2.5 m=H/2; ξ=1, (c) θv=1.0 m=H/5; ξ=6, (d) θv=1.0 m=H/5; ξ=12.
Results and evaluation
Figure 5 compares deterministic finite element and analytical solutions for factor of safety (FOS), for three distributions of cu: cu constant with depth (Taylor 1937); cu increasing linearly with depth (Hunter & Schuster 1968); and cu increasing linearly, with depth, from zero at z=0 m (Gibson & Morgenstern 1962). In each case, the FOS has been found for a slope with an average cu of 40 kPa, i.e. over the 5 m depth. Hence, the three cu distributions are: 40 kPa, 20–60 kPa and 0–80 kPa. Note that, for the latter two cases, cu increases linearly relative to z=0 m, i.e. the horizontal ground surface (not the slope face). Hence, the analyses are applicable to cuttings in various types of clay.
Typical random fields of cu for μ∝depth. (a) θv=1.0 m=H/5; ξ=1, (b) θv=1.0 m=H/5; ξ=6.
Comparing deterministic finite element and analytical solutions for factors of safety.
Each computed curve has been produced by conducting a series of analyses. Hence, for a given curve, each data point represents a separate analysis, in which the crest settlement (Δ) has been computed by applying gravitational forces in a single increment. Starting from FOS=1.0, the cu values have been decreased from one analysis to the next, by dividing the original cu distribution by progressively larger values of FOS, until the analysis fails to converge in a reasonable number of equilibrium iterations. This gives the factor of safety corresponding to an average cu of 40 kPa.
Figure 5 shows good agreement with theory. As expected, for the same average undrained shear strength, a greater factor of safety is obtained when cu increases with depth, due to failure being dependent on slope height. However, note that for the special case of cu increasing from zero, slope stability depends only on slope angle: i.e. it is independent of H, and the slope is liable to fail en-masse, rather than along a well-defined rupture surface originating from the slope toe. Figure 5 indicates that the specified iteration limit of 500 has been suitable for detecting collapse in each of the three cases.
Figure 6 shows typical stochastic responses for mean cu constant with depth: these are based on μ=40 kPa, as in Figure 5, and COV=0.3. The results show that the mean stochastic FOS is lower than the deterministic solution, and that the standard deviation of the output statistics increases with both increasing θ and increasing degree of anisotropy.
Influence of θ on stochastic solutions for μ constant with depth. (a) θv=0.25 m=H/20; ξ=1, (b) θv=1.0 m=H/5; ξ=1, (c) θv=1.0 m=H/5; ξ=6.
Paice & Griffiths (1997) suggested a simple procedure for estimating slope reliability, R. That is, for a given set of statistics, multiple realizations are performed, in which each realization involves the application of gravity in a single increment. R is then the percentage of realizations which do not reach failure, i.e.
in which N is the total number of realizations in the ensemble, and Nf is the number of realizations reaching failure (as indicated by non-convergence in 500 iterations). This technique has been used for obtaining the results shown in Figures 7 and 8 , in which 300realizations have been carried out for each set of statistics.
Influence of θ on reliability for various factors of safety and ξ=1. (a) μ constant with depth, (b) μ∝depth.
Firstly, Figure 7 shows the influence of θ on reliability, for isotropic heterogeneity and various factors of safety. In contrast to the method adopted in Figures 5 and 6, FOS=1.0 corresponds to a slope at the point of failure, based on deterministic analysis. Hence, the mean cu for limit equilibrium has been found by dividing the cu distribution used in Figure 5 by the FOS computed in the same figure: i.e. in Figure 7(a), for mean cu (and COV) constant with depth (cf. Taylor 1937), μ=40.0/2.46=16.3 kPa. Results for higher factors of safety have then been obtained by factoring μ (=16.3 kPa) and σ (=4.9 kPa) by FOS, i.e. when computing distributions of cu using equation (6). The same procedure has been adopted in Figure 7(b), for the case of the mean and standard deviation of cu increasing linearly from zero at z=0 m, and COV=0.3 (cf. Gibson & Morgenstern 1962).
Figure 7(a) shows that, for mean cu constant with depth, FOS≥1.4 is needed for there to be virtually no chance of failure for θ/H≤0.2, while, for mean cu proportional to depth, a higher FOS (≥1.8) is required. This is due to the greater range of possible rupture surfaces, as shown by the shear strain invariant contours plotted in Figures 9 and 10 (note that lighter colours indicate larger strains). For constant μ (Fig. 9), all rupture surfaces originate from the slope toe (or very near to it), due to the slope height having a greater influence on stability than variations in cu. Therefore, heterogeneity mainly influences the shape of a rupture surface and the location at which it exits the ground surface. In contrast, for μ∝z and COV constant with depth (Fig. 10), heterogeneity means that failure can originate from any point along the slope, due to stability being independent of H for this case.
Influence of ξ on reliability for various factors of safety and θv/H=0.2. (a) μ constant with depth, (b) μ∝depth.
Typical failure modes for μ constant with depth. (a) θv=1.0 m=H/5; ξ=1, (b) θv=1.0 m=H/5; ξ=3, (c)– (d) θv=1.0 m=H/5; ξ=6 (results of 2 realizations).
Typical failure modes for μ∝depth. (a)– (c) θv=1.0 m=H/5; ξ=1 (results of 3 realizations), (d) θv=1.0 m=H/5; ξ=12.
Figures 7(a) and 7(b) show that, when FOS=1.0, R is lower for small values of θ. This is most easily explained with reference to Figures 6(a) and 6(b). For small values of θ, the standard deviation of the output statistics is much smaller than for larger values of θ. Therefore, as the mean response is lower than the deterministic solution (for all values of θ), there is a decreased probability of a stable solution when θ is small. However, it takes a relatively small increase in FOS to push the mean response above the deterministic solution. When this happens, R increases rapidly for small values of θ, as seen in Figure 7, due to the small standard deviation now making the possibility of failure unlikely. Hence, for intermediate and higher factors of safety, reliability decreases as θ increases.
Figure 8 shows the influence of horizontal anisotropy (of the heterogeneity) on reliability, for θv/H=0.2. For mean cu constant with depth (Fig. 8(a)), R increases with increasing degree of anisotropy for lower factors of safety, due to the larger range of possible solutions (compare Figures 6(b) and 6(c)). However, the rate at which R increases with increasing FOS is slower than for isotropic heterogeneity, due to the larger standard deviation of the output statistics. Indeed, Figure 8(a) shows that, for ξ>1.0, a larger FOS is needed for R≈100%. A slower increase in reliability with FOS is also observed for μ∝z (Fig. 8(b)). In this case, for FOS=1.0, there is a negligible increase in reliability with increasing ξ.
Conclusions
Stochastic analysis leads to reliability, i.e. the probability that failure will not occur; a more meaningful definition of stability than factor of safety, for which there is no information regarding probability of failure. However, stochastic analysis does have a role to play in conventional (deterministic) design, as factors of safety corresponding to certain levels of reliability (e.g. 95%) may be derived, or, alternatively, characteristic values determined. In this case, spatial variability, as well as pointwize variability, should be considered: cf. Eurocode 7 (DD ENV 1997-1 1995), which suggests characteristic values based only on pointwize variability, due largely to insufficient knowledge in quantifying and analysing heterogeneity.
This paper has applied stochastics to a classical slope stability problem, using typical statistics found in literature. For low and intermediate values of the coefficient of variation of cu (e.g. 0.1–0.3), as reported by many researchers, a normal pdf may be used to describe pointwize variability. Data for spatial variability are scarce. However, evidence suggests that the vertical scale of fluctuation will often be small, relative to the horizontal scale of fluctuation and to the height of a slope. For a given factor of safety, reliability decreases for problems in which there are a greater range of possible rupture surfaces, as is the case for mean cu proportional to depth.
Anisotropy (of the heterogeneity) can either promote or hinder slope stability. For the horizontal anisotropy considered in this paper, modest changes in reliability have been determined. However, Hicks & Onisiphorou (2000) and Onisiphorou (2000) have demonstrated that anisotropy can sometimes have a significant and destabilizing influence: in particular, for problems in which anisotropy is inclined at some angle to the horizontal, as might be expected in natural slopes and man-made fills; also, for materials with a potential for dramatic strength loss, e.g. loose hydraulic fills and sensitive clays. In such cases, destabilization may result from failure mechanisms propagating along semi-continuous weak zones.
- © 2002 The Geological Society of London
