# Sampling Efficiency in Spatially Varying Soils for Slope Stability Assessment

Authors: Yajun Li; Cheng Qian; Kang Liu
Journal: Advances in Civil Engineering
 (2019)

---

## Abstract

Site investigations are usually carried out in geotechnical engineering to track the range of design parameters. Due to the inherent soil spatial variability and usually limited scope of site exploration programs, the design parameters are usually uncertain at locations where test samples are not taken. This uncertainty often propagates to the response of geotechnical structures such as soil slopes. This paper developed a conditional simulation framework to investigate the sampling efficiency (i.e., sampling location and sampling distance) in designing slopes in spatially varying soils. A performance-based sampling efficiency index is proposed to achieve this. It is found that the optimal location to take vertical samples are in the vicinity of slope crest to slope midpoint and the optimal distance within a limited exploration scope is such that the additional sample subdivides the influence domain to either side of the first vertical sample.


---

## Body

## 1. Introduction

Spatial variability (i.e., heterogeneity) of soil properties can have significant effects on geotechnical structures [1]. Explicit consideration of spatial variability often optimises the design of such structures [2]. Therefore, it is important to carry out a site investigation with a certain scope (maybe in multiple stages) and collect and process the measurement data statistically to characterise the variability [3–13]. The statistical outcome, e.g., the distribution type parametrized by the mean soil property value and the coefficient of variation (i.e., the ratio of standard deviation to the mean), and the spatial correlation structure, may be used as inputs to a random field model capable of simulating the spatial variation.Geotechnical investigations include field tests and laboratory experiments. For example, cone penetration test (CPT) is frequently performed in situ and triaxial tests carried out in the laboratory. In view of estimating the spatial correlation structure, CPT measurements are often preferable to conventional laboratory tests as the amount of data is usually much larger. For example, a database made of CPT measurements from Oslo, Norway, was used to estimate the vertical correlation statistics [5], and a database comprising CPT measurements from Adelaide, Australia, was used to estimate both the vertical and horizontal correlation distances. Thus, this paper will focus on CPT sampling.An obvious question arises regarding the design of the sampling strategy (i.e., sampling location and sampling distance between measurements) in the soil deposit for a particular geostructure. For example, if an embankment slope is to be constructed in a soil, a number of CPTs are proposed to be conducted along a straight line in the longitudinal direction (e.g.,x); where is the best location to take the CPTs in the perpendicular (e.g., y−z) plane such that the performance of the proposed slope would have the smallest uncertainty, e.g., in terms of stability uncertainty indicated by the standard deviation of realised factor of safety? This paper is devoted to answering this question.Li et al. [2] investigated the optimal number of CPTs in the longitudinal direction with regard to slope stability, without addressing phased site investigations in that direction. For instance, after a first phase site test with one CPT in the middle, what would be the best location to carry out the second phase investigation if it is decided to carry out a second phase test? This question is not answered there. Although in 2D, the current study complements the previous study by looking at the sampling distance in the cross section (i.e., perpendicular to the longitudinal direction), in a more often multistaged site investigation. This study focuses on a two-staged investigation; as in slope stability applications of this kind, the scope in the y−z plane is quite often limited to no more than 3 CPT profiles (more CPTs are often taken in the third dimension). In particular, this study generalises the previous findings by investigating more site spatial variability scenarios (i.e., with different horizontal variation situations) in the cross section, although it assumes an infinite scale of fluctuation in the third dimension (which implies an infinite number of identical CPTs in the third dimension). In the previous research, when trying to find the best location for the second row of CPTs in the cross section, the best sampling interval relative to the first row of CPTs happens to be half of the correlation distance considered. By investigating more spatial variability scenarios, these conclusions will be generalised in this paper to avoid misunderstanding of the previous findings.For reducing the uncertainty in slope stability assessment, Yang et al. [14] used a framework based on a hypothesis test to investigate the probability of Type I and Type II errors in slope stability evaluation (i.e., the first type of error is that the slope is actually unsafe, but the slope stability analysis suggests that the slope is safe. This is termed Type I error. The second type of error is that the slope is actually safe, but the slope stability analysis suggests that the slope is unsafe. This is termed Type II error). Similar to [2], they found the optimal location to carry out CPT site investigations to be roughly at the slope crest, although they used a different method of analysis. However, they did not look at the other important aspect of efficient sampling, i.e., the optimal distance between samples. Besides, they used kriging to estimate the unknown property field, and this was then plugged into a finite element model to see if the slope with the kriged property is stable or not, compared to a slope with a reference random field property. However, it is important to understand the difference between estimation (e.g., the kriged field) and random field simulation (e.g., the random field). Kriging tends to smooth out details and extreme values of the original data set at places other than the measurement data locations. This means that small values are overestimated, while large values tend to be underestimated. Thus, the kriging estimator does not reproduce the true local variability of the property, whereas conditional random field simulation does. Therefore, this paper implemented a conditional random field simulator, combining the unconditional random field and kriging estimation, to investigate both aspects (i.e., optimal location and optimal distances) of efficient sampling with regard to slope stability.With application to two clay slope stability examples, it is demonstrated that the best locations for carrying out CPT tests can be identified by the proposed approach, thereby to increase confidence in a slope’s stability state (stable or not). The first example seeks to find the optimum locations for site investigations in the first phase, before moving on to find the optimal distance between tests for various site variabilities. The second example considers a slope with a foundation layer to see if the inclusion of a soil foundation layer changes the general findings. The idea is to further generalise the findings regarding the optimal locations and optimal distances of tests in the first example.For simplicity, this paper focuses on applications involving only a single soil layer (i.e., a single layer characterised by a statistically homogeneous undrained shear strength), although the extension to multiple soil layers is straightforward. Moreover, the effect of random variation in the boundary locations between different soil layers can also be easily incorporated by conditioning to known boundary locations (e.g., corresponding to where the CPTs have been carried out). For example, Li et al. [15] considered the influence of sampling layouts on slope reliability using a coupled Markov chain (CMC) model (i.e., a sort of conditional simulation model [16]), Deng et al. [17] then extended the framework to take account of soil spatial variability within the different soil geological units. However, the test borehole data have not been fully utilised with regards to the spatial distribution of the boreholes (when simulating the soil spatial variability within each soil unit). Also, one important and related issue, that is, the relationship between sampling interval and fluctuation scale of soil geological units, was not addressed in their work. The framework proposed in this paper can be extended to include soil geological stratigraphy uncertainty and is a subject of ongoing research.


## 2. Conditional Random Fields

This study uses the method of Li et al., and the method description partly reproduces their wording [2]. A conditional random field, constrained by the known measurement values at some locations, can be obtained from three different fields [18–20]:(1)Zrcx=Zrux+Zkmx−Zksx,where x denotes the spatial location, Zrcx is a random field conditional upon the measurement, Zrux is an unconditionally simulated random field, Zkmx is the kriged field based on measurement values at xi (i=1,2,…,N, where N is the number of measurement locations), and Zksx is the kriged field based on the randomly simulated values at the same positions xi (i=1,2,…,N) as the measurements.Several methods can be used to simulate the unconditional random field [21–23]. The local average subdivision (LAS) method [24] is used in this paper. The local average subdivision process [24, 25] proceeds in a top-down recursive manner. LAS in 2D involves the subdivision of a square parent cell into 4 equally sized smaller cells. Initially, a single square cell is used to represent the random field domain. It is then uniformly subdivided into smaller square cells. This goes on stage by stage (i.e., a stage means a subdivision level) until the required level of subdivision. An existing LAS program implemented by Spencer [26] and later improved by Li et al. [1, 27, 28] has been used in this paper. The program uses the covariance function as follows:(2)Cτ1,τ2=σ2exp−2τ1θ1−2τ2θ2,where θ1 and θ2 denote the scale of fluctuation (SOF) in the lateral and vertical coordinate directions, respectively, τ1 and τ2 are the lag distances in the two respective directions, and σ2 is the soil property variance. However, due to the subdivision algorithm itself not being capable of preserving correlation anisotropy [21], an isotropic random field is initially generated in this implementation, and this field is then squashed and/or stretched in the vertical and/or horizontal directions, respectively, to give the target anisotropic random field.Kriging [29], which is commonly used in geostatistics [30, 31], has been used to obtain the kriged fields. The following section is devoted to a brief review of the method and its underlying theory.


### 2.1. Kriging Interpolation Theory

Kriging interpolation can be used to obtain a best linear unbiased prediction of soil properties (Z) by using known measurement data [31, 32] and estimated correlation structure and assuming a stationary mean and spatial covariances (or variograms). Suppose that Z1,Z2,…,ZN are observations of the random field Zx at points x1,x2,…,xN (i.e., Zi=Zxi,i=1,2,…,N). The best linear unbiased estimate Z^x0 of the soil property at some target location x0 is given by(3)Z^x0=∑i=1NλiZi=∑i=1Nλix0Zxi,where N is the total number of observations and λi denotes the unknown weighting factor (to be determined) associated with observation point xi. The kriging estimate is a weighted average of the values of Zi. The weights are required to sum up to 1 to ensure an unbiased estimation, i.e.,(4)∑i=1Nλi=1.The weightsλi in equation (3), for the estimation at any point x0, can be solved by finding the minimum of the variance of the kriging error σe2.The Lagrange method can be used [31] to minimise the error variance. This leads to a system of N+1 equations with unknowns λi and μ, for a constant mean:(5)γx1−x1⋯γx1−xN1⋮⋱⋮⋮γxN−x1⋯γxN−xN11⋯10λ1⋮λNμ=γx1−x0⋮γxN−x01,where a Lagrangian parameter μ is introduced and γτ is the variogram. It is related to the covariance function Cτ (equation (2)) by(6)γτ=C0−Cτ=σ2−Cτ.The error variance is then expressed as(7)σe2=μ+∑i=1Nλiγxi−x0−γ0,where γ0=0 according to equation (6).For a mean following some trend, the modification to equation (5) is straightforward, and interested readers are referred to Chapter 4 in [32].Equation (5) may be expressed in compact form as(8)γlhsλx=γrhsand solved by inverting the left-hand side matrix, i.e.,(9)λx=γlhs−1γrhs.Then, the estimation error variance (according to equation (7)) can be expressed by(10)σe2=γlhs−1γrhsTγrhs,where λi is essentially a function of the configuration (i.e., relative positioning) of measurement points xi and the target point x0.Note that the left-hand side matrixγlhs is only a function of covariances or variograms between pairs of the observation points. Therefore, for a certain configuration of measurement locations, it suffices to invert once, and then equations (9) and (3) can be used repeatedly for building up the field of interest (i.e., the kriged field), i.e., through point by point estimation in space. In contrast, as the target spatial point x0 moves from point to point in space, the right-hand side vector γrhs changes accordingly. Thus, different weight vectors λx are obtained and used in equation (3) to estimate point by point in order to form the kriged field.


### 2.2. Kriging Implementation

In geotechnical engineering, it is common for a sampling strategy to follow some pattern [3]. For example, a regular grid of CPT sampling on the ground surface is often planned [7]. Although kriging can be easily implemented based on data point coordinates (whether they follow a regular pattern or not), it generally requires an input file containing all the data point coordinates. In the case of some sampling design with a regular pattern, the required input can, however, be simplified and the implementation to compensate this input simplification can be nontrivial. The most fundamental part of the algorithm is to form the left-hand side matrix in equation (8). This paper shows a possible implementation of kriging in the case when the sampling strategy follows a regular grid as is shown in Figure 1. A global numbering scheme is first set up, the left-hand side matrix is then constructed, and the right-hand side vector and the unknown weighting vector are formed.Figure 1
Example illustration of CPT sampling strategy (k = 2, m = 5).Letk×m CPTs be located in the form of a rectangular grid (plan view, x−y). That is, k CPTs are lined up in the x direction and m CPT profiles in the y direction (Figure 1). Assuming that for each CPT profile, there are n data points. In order to set up the relevant matrices and vectors, a global point or cell numbering scheme is proposed for all the CPT measurement points first, and this is shown in Figure 2 for the case of k = 1 (i.e., when looking at a cross section y−z). Although the investigation on slope stability design later in the paper will be mainly focusing on 2D analysis, the following implementation is formulated for the general case of 3D CPT arrangement. For the special case of cross-sectional analysis, k=1.Figure 2
Example CPT data grid (k = 1): global numbering of data points at a cross section.


#### 2.2.1. Forming LHS Matrix

Following the basic equation (equation (8)) of size N+1=kmn+1, the left-hand side matrix is constructed as(11)γlhs=v1,1v1,2v1,3…v1,mv1,m+1v1,m+2v1,m+3…v1,2m…v1,k−1m+1v1,k−1m+2v1,k−1m+3…v1,km1v2,1v2,2v2,3…v2,mv2,m+1v2,m+2v2,m+3…v2,2m…v2,k−1m+1v2,k−1m+2v2,k−1m+3…v2,km1v3,1v3,2v3,3…v3,mv3,m+1v3,m+2v3,m+3…v3,2m…v3,k−1m+1v3,k−1m+2v3,k−1m+3…v3,km1⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋱⋮⋮⋮⋱⋮⋮vkm,1vkm,2vkm,3…vkm,mvkm,m+1vkm,m+2vkm,m+3…vkm,2m…vkm,k−1m+1vkm,k−1m+2vkm,k−1m+3…vkm,km1111…1111…1…111…10,where vi,j is a submatrix that represents the correlation structure between CPTi and CPTj (where each CPT has n data points):(12)vi,j=di−1n+1,j−1n+1di−1n+1,j−1n+2…di−1n+1,j−1n+ndi−1n+2,j−1n+1di−1n+2,j−1n+2…di−1n+2,j−1n+ndi−1n+3,j−1n+1di−1n+3,j−1n+2…di−1n+3,j−1n+n⋮⋮⋱⋮di−1n+n,j−1n+1di−1n+n,j−1n+2…di−1n+n,j−1n+n,where i,j=1,2,3,…,m,m+1,m+2,m+3,…,2m,…,k−1m+1,k−1m+2,k−1m+3,…,km and dr,s are the components of the submatrix vi,j (where r=i−1n+1,…,i−1n+n; s=j−1n+1,…,j−1n+n) and can be expressed in the form of variograms between data points r and s (equations (6) and (2)).


#### 2.2.2. Forming RHS Vector and the Weighting Vector

The right-hand side vector is organised as(13)γrhs=v1v2v3…vkm1T,where vp is a subvector representing the correlation structure between the target estimation point and CPTp:(14)vp=dp−1n+1dp−1n+2dp−1n+3…dp−1n+nT,where p=1,2,3,…,m,m+1,m+2,m+3,…,2m,…,k−1m+1,k−1m+2,k−1m+3,…,km and dt (where t=p−1n+1,…,p−1n+n) are the components of the subvector vp, which can also be expressed in the form of variograms (equation (6) and (2)) between data points t and the target point (Figure 2).The unknown weight vector is formed as(15)λx=λ1λ2λ3…λkmμT,where λq is the weight subvector for CPTq:(16)λq=λq−1n+1λq−1n+2λq−1n+3…λq−1n+nT,where q=1,2,3,…,m,m+1,m+2,m+3,…,2m,…,k−1m+1,k−1m+2,k−1m+3,…,km.Note that equation (13) is formed and the unknown vector is solved cell by cell (Figure 2). The blank cells in Figure 2 denote all the cell values that are to be estimated (i.e., one by one), whereas the grey cells denote the CPT measurements. Thus, the series of discrete values estimated are assigned to the blank square cells to represent the best linear estimation of the soil property at those locations.


## 2.1. Kriging Interpolation Theory

Kriging interpolation can be used to obtain a best linear unbiased prediction of soil properties (Z) by using known measurement data [31, 32] and estimated correlation structure and assuming a stationary mean and spatial covariances (or variograms). Suppose that Z1,Z2,…,ZN are observations of the random field Zx at points x1,x2,…,xN (i.e., Zi=Zxi,i=1,2,…,N). The best linear unbiased estimate Z^x0 of the soil property at some target location x0 is given by(3)Z^x0=∑i=1NλiZi=∑i=1Nλix0Zxi,where N is the total number of observations and λi denotes the unknown weighting factor (to be determined) associated with observation point xi. The kriging estimate is a weighted average of the values of Zi. The weights are required to sum up to 1 to ensure an unbiased estimation, i.e.,(4)∑i=1Nλi=1.The weightsλi in equation (3), for the estimation at any point x0, can be solved by finding the minimum of the variance of the kriging error σe2.The Lagrange method can be used [31] to minimise the error variance. This leads to a system of N+1 equations with unknowns λi and μ, for a constant mean:(5)γx1−x1⋯γx1−xN1⋮⋱⋮⋮γxN−x1⋯γxN−xN11⋯10λ1⋮λNμ=γx1−x0⋮γxN−x01,where a Lagrangian parameter μ is introduced and γτ is the variogram. It is related to the covariance function Cτ (equation (2)) by(6)γτ=C0−Cτ=σ2−Cτ.The error variance is then expressed as(7)σe2=μ+∑i=1Nλiγxi−x0−γ0,where γ0=0 according to equation (6).For a mean following some trend, the modification to equation (5) is straightforward, and interested readers are referred to Chapter 4 in [32].Equation (5) may be expressed in compact form as(8)γlhsλx=γrhsand solved by inverting the left-hand side matrix, i.e.,(9)λx=γlhs−1γrhs.Then, the estimation error variance (according to equation (7)) can be expressed by(10)σe2=γlhs−1γrhsTγrhs,where λi is essentially a function of the configuration (i.e., relative positioning) of measurement points xi and the target point x0.Note that the left-hand side matrixγlhs is only a function of covariances or variograms between pairs of the observation points. Therefore, for a certain configuration of measurement locations, it suffices to invert once, and then equations (9) and (3) can be used repeatedly for building up the field of interest (i.e., the kriged field), i.e., through point by point estimation in space. In contrast, as the target spatial point x0 moves from point to point in space, the right-hand side vector γrhs changes accordingly. Thus, different weight vectors λx are obtained and used in equation (3) to estimate point by point in order to form the kriged field.


## 2.2. Kriging Implementation

In geotechnical engineering, it is common for a sampling strategy to follow some pattern [3]. For example, a regular grid of CPT sampling on the ground surface is often planned [7]. Although kriging can be easily implemented based on data point coordinates (whether they follow a regular pattern or not), it generally requires an input file containing all the data point coordinates. In the case of some sampling design with a regular pattern, the required input can, however, be simplified and the implementation to compensate this input simplification can be nontrivial. The most fundamental part of the algorithm is to form the left-hand side matrix in equation (8). This paper shows a possible implementation of kriging in the case when the sampling strategy follows a regular grid as is shown in Figure 1. A global numbering scheme is first set up, the left-hand side matrix is then constructed, and the right-hand side vector and the unknown weighting vector are formed.Figure 1
Example illustration of CPT sampling strategy (k = 2, m = 5).Letk×m CPTs be located in the form of a rectangular grid (plan view, x−y). That is, k CPTs are lined up in the x direction and m CPT profiles in the y direction (Figure 1). Assuming that for each CPT profile, there are n data points. In order to set up the relevant matrices and vectors, a global point or cell numbering scheme is proposed for all the CPT measurement points first, and this is shown in Figure 2 for the case of k = 1 (i.e., when looking at a cross section y−z). Although the investigation on slope stability design later in the paper will be mainly focusing on 2D analysis, the following implementation is formulated for the general case of 3D CPT arrangement. For the special case of cross-sectional analysis, k=1.Figure 2
Example CPT data grid (k = 1): global numbering of data points at a cross section.


### 2.2.1. Forming LHS Matrix

Following the basic equation (equation (8)) of size N+1=kmn+1, the left-hand side matrix is constructed as(11)γlhs=v1,1v1,2v1,3…v1,mv1,m+1v1,m+2v1,m+3…v1,2m…v1,k−1m+1v1,k−1m+2v1,k−1m+3…v1,km1v2,1v2,2v2,3…v2,mv2,m+1v2,m+2v2,m+3…v2,2m…v2,k−1m+1v2,k−1m+2v2,k−1m+3…v2,km1v3,1v3,2v3,3…v3,mv3,m+1v3,m+2v3,m+3…v3,2m…v3,k−1m+1v3,k−1m+2v3,k−1m+3…v3,km1⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋱⋮⋮⋮⋱⋮⋮vkm,1vkm,2vkm,3…vkm,mvkm,m+1vkm,m+2vkm,m+3…vkm,2m…vkm,k−1m+1vkm,k−1m+2vkm,k−1m+3…vkm,km1111…1111…1…111…10,where vi,j is a submatrix that represents the correlation structure between CPTi and CPTj (where each CPT has n data points):(12)vi,j=di−1n+1,j−1n+1di−1n+1,j−1n+2…di−1n+1,j−1n+ndi−1n+2,j−1n+1di−1n+2,j−1n+2…di−1n+2,j−1n+ndi−1n+3,j−1n+1di−1n+3,j−1n+2…di−1n+3,j−1n+n⋮⋮⋱⋮di−1n+n,j−1n+1di−1n+n,j−1n+2…di−1n+n,j−1n+n,where i,j=1,2,3,…,m,m+1,m+2,m+3,…,2m,…,k−1m+1,k−1m+2,k−1m+3,…,km and dr,s are the components of the submatrix vi,j (where r=i−1n+1,…,i−1n+n; s=j−1n+1,…,j−1n+n) and can be expressed in the form of variograms between data points r and s (equations (6) and (2)).


### 2.2.2. Forming RHS Vector and the Weighting Vector

The right-hand side vector is organised as(13)γrhs=v1v2v3…vkm1T,where vp is a subvector representing the correlation structure between the target estimation point and CPTp:(14)vp=dp−1n+1dp−1n+2dp−1n+3…dp−1n+nT,where p=1,2,3,…,m,m+1,m+2,m+3,…,2m,…,k−1m+1,k−1m+2,k−1m+3,…,km and dt (where t=p−1n+1,…,p−1n+n) are the components of the subvector vp, which can also be expressed in the form of variograms (equation (6) and (2)) between data points t and the target point (Figure 2).The unknown weight vector is formed as(15)λx=λ1λ2λ3…λkmμT,where λq is the weight subvector for CPTq:(16)λq=λq−1n+1λq−1n+2λq−1n+3…λq−1n+nT,where q=1,2,3,…,m,m+1,m+2,m+3,…,2m,…,k−1m+1,k−1m+2,k−1m+3,…,km.Note that equation (13) is formed and the unknown vector is solved cell by cell (Figure 2). The blank cells in Figure 2 denote all the cell values that are to be estimated (i.e., one by one), whereas the grey cells denote the CPT measurements. Thus, the series of discrete values estimated are assigned to the blank square cells to represent the best linear estimation of the soil property at those locations.


## 2.2.1. Forming LHS Matrix

Following the basic equation (equation (8)) of size N+1=kmn+1, the left-hand side matrix is constructed as(11)γlhs=v1,1v1,2v1,3…v1,mv1,m+1v1,m+2v1,m+3…v1,2m…v1,k−1m+1v1,k−1m+2v1,k−1m+3…v1,km1v2,1v2,2v2,3…v2,mv2,m+1v2,m+2v2,m+3…v2,2m…v2,k−1m+1v2,k−1m+2v2,k−1m+3…v2,km1v3,1v3,2v3,3…v3,mv3,m+1v3,m+2v3,m+3…v3,2m…v3,k−1m+1v3,k−1m+2v3,k−1m+3…v3,km1⋮⋮⋮⋱⋮⋮⋮⋮⋱⋮⋱⋮⋮⋮⋱⋮⋮vkm,1vkm,2vkm,3…vkm,mvkm,m+1vkm,m+2vkm,m+3…vkm,2m…vkm,k−1m+1vkm,k−1m+2vkm,k−1m+3…vkm,km1111…1111…1…111…10,where vi,j is a submatrix that represents the correlation structure between CPTi and CPTj (where each CPT has n data points):(12)vi,j=di−1n+1,j−1n+1di−1n+1,j−1n+2…di−1n+1,j−1n+ndi−1n+2,j−1n+1di−1n+2,j−1n+2…di−1n+2,j−1n+ndi−1n+3,j−1n+1di−1n+3,j−1n+2…di−1n+3,j−1n+n⋮⋮⋱⋮di−1n+n,j−1n+1di−1n+n,j−1n+2…di−1n+n,j−1n+n,where i,j=1,2,3,…,m,m+1,m+2,m+3,…,2m,…,k−1m+1,k−1m+2,k−1m+3,…,km and dr,s are the components of the submatrix vi,j (where r=i−1n+1,…,i−1n+n; s=j−1n+1,…,j−1n+n) and can be expressed in the form of variograms between data points r and s (equations (6) and (2)).


## 2.2.2. Forming RHS Vector and the Weighting Vector

The right-hand side vector is organised as(13)γrhs=v1v2v3…vkm1T,where vp is a subvector representing the correlation structure between the target estimation point and CPTp:(14)vp=dp−1n+1dp−1n+2dp−1n+3…dp−1n+nT,where p=1,2,3,…,m,m+1,m+2,m+3,…,2m,…,k−1m+1,k−1m+2,k−1m+3,…,km and dt (where t=p−1n+1,…,p−1n+n) are the components of the subvector vp, which can also be expressed in the form of variograms (equation (6) and (2)) between data points t and the target point (Figure 2).The unknown weight vector is formed as(15)λx=λ1λ2λ3…λkmμT,where λq is the weight subvector for CPTq:(16)λq=λq−1n+1λq−1n+2λq−1n+3…λq−1n+nT,where q=1,2,3,…,m,m+1,m+2,m+3,…,2m,…,k−1m+1,k−1m+2,k−1m+3,…,km.Note that equation (13) is formed and the unknown vector is solved cell by cell (Figure 2). The blank cells in Figure 2 denote all the cell values that are to be estimated (i.e., one by one), whereas the grey cells denote the CPT measurements. Thus, the series of discrete values estimated are assigned to the blank square cells to represent the best linear estimation of the soil property at those locations.


## 3. Conditional Simulation

The conditional simulation of geotechnical performance based on finite elements and the above kriging-based conditional random fields is presented here to show the workflow of the procedure. A flowchart for undertaking such a simulation is shown in Figure3. The procedure first involves kriging a soil property field based on real measurement data using the previous implementation. This step only needs to be done once. It then proceeds to generate the conditional random field realisations of the soil property, which is, as described in the previous section, comprised of the unconditional random fields and kriged fields based on the randomly simulated values at those same locations as those of the real measurements. The third step goes on with a finite element analysis carried out repeatedly for the desired number of Monte Carlo realisations, by taking the generated conditional random field values at the Gaussian point level for each realisation. Unlike the first step, the second and third steps fall into the realisation loop to simulate the possible realisations of conditional random fields and their corresponding structural performance in a finite element analysis. Note that the procedure is herein generally presented for various possible geotechnical applications, although the following sections are devoted to the particular problem of slope stability assessment. The procedure is equally applicable to other ultimate limit state design problems as well as serviceability limit state design problems where, for instance, the displacement response may be of concern [33].Figure 3
The procedure for conditional RFEM simulation (#1: via equations (5) and (3); #2: via equation (1)).Note that in order to use the direct CPT data, i.e., cone resistance and sleeve friction, in the analyses reported here, a conversion or transformation model is needed to relate the test measurement (e.g., cone resistance) to an appropriate design property (e.g., the undrained shear strength) [34]. For simplicity, this paper only shows how such converted or transformed data can be used. However, the framework presented here can be easily extended to include measurement and transformation errors. Should those additional uncertainties be included in the model, the estimation variance at sampled locations would have been larger than zero, although they would have been still smaller than those at locations that are further away from samples.To show the preservation of the measured values within each realisation, and to demonstrate the effect of conditioning the random fields, samples are to be taken in a 5 m high (z) and 10 m wide (y) clay block, of which the undrained shear strength cu is spatially varying. The block is discretised into 40 × 80 square cells, each of side dimension 0.125 m. The mean of cu is 20 kPa, and the standard deviation of cu is 6 kPa. The scale of fluctuation of the strength heterogeneity is θh = 2.0 m and θv = 1.0 m. The coefficient of variation of 0.3 and the chosen values of the scale of fluctuation are well within the recommended ranges for undrained shear strength [35]. CPT measurements are taken in the y direction at two locations, with each CPT comprising n = 40 data points at 0.125 m vertical interval. Note that a single independent realisation of the spatially varying strength (i.e., representing the “actual” or “real” in situ variability) is artificially generated first and sampled in a vertical line to obtain these “measured” CPT data. The first CPT profile is at y = 2.5 m, and the horizontal spacing between CPTs is Δy = 3.75 m.An example realisationi is shown in Figure 4, to illustrate the various fields involved in generating the conditional random field. It shows (a) the conditional random field, (b) the kriged field based on the measured data (taken from the reference field), (c) the unconditional random field generated using LAS, and (d) the kriged field based on the randomly simulated cell values (from field (c)) at the exact measurement locations. The location from which the CPTs are taken is also shown in Figure 4, together with the measured (known) CPT profiles in Figures 4(a) and 4(b). It is seen that the conditional field honors the measured CPT data at the measurement locations, eliminating the unrealistic values (larger or smaller) from the unconditional fields (e.g., those values that correspond to the area of the circle shown in the figure). Note that, darker and lighter colours represent low and high values, respectively, and in order to better compare the various fields, a global colour scale is used for all part figures in Figures 4 and 5.Figure 4
The four fields for realisationi. (a) Conditional random field i. (b) Kriged field based on measurements. (c) Unconditional random field i. (d) Kriged field i based on random values at measurement locations.
(a)
(b)
(c)
(d)Figure 5
The four fields for realisationj. (a) Conditional random field j. (b) Kriged field based on measurements. (c) Unconditional random field j. (d) Kriged field j based on random values at measurement locations.
(a)
(b)
(c)
(d)Figure5 shows another example realisation j, with the same series of part figures as in Figure 4, to show the differences between realisations. The difference between Figures 4(c)–4(d) and 5(c)–5(d) with regard to relative locations of weak and strong soil pockets is apparent, causing the difference between Figures 4(a) and 5(a), albeit less apparent due to the conditioning upon the two CPT profiles. Note that Figure 5(b) is simply a replication of Figure 4(b) to facilitate visualisation and comparison, as mentioned previously, this field only needs to be kriged once outside the realisation loop in contrast to the kriged field in Figures 4(d) and 5(d).


## 4. Sampling Efficiency Applied to Slope Cut Design

In order to quantify and assess the efficiency of a sampling design for slope stability, a performance-based sampling efficiency index, based on the standard deviations of conditional and unconditional simulation, is here defined as(17)Ise=σuσi,where σu and σi are the standard deviation of the realised factor of safety for the unconditional and conditional slope stability simulation when the sampling location is at position i. Hence, Ise=1 for an unconditioned simulation. Note that Ise may be larger or smaller than 1.0, depending the kriging error variance (see equation (10)) which could be larger than the input variance of soil property due to ordinary kriging also estimating the unknown mean (in this case, the conditional variance for points that are far away from the measurement points would be larger than the unconditional variance).Two simple slope stability examples are presented in this section, in order to demonstrate the capability of conditional simulation as a useful aid to geotechnical sampling design. In particular, the optimum locations and sampling distances for CPT profiles are investigated in assessing the stability of the first slope, in order to minimise the stability uncertainty of the slope founded on a firm base. The same sampling procedure is investigated in the second example, involving a slope with a foundation layer to investigate the influence of including a foundation layer.In the two examples, results are both presented in terms of the uncertainty in the slope stability (with respect to the realised factor of safety). The strength reduction method [36] in finite elements has been used to obtain the factors of safety. The analyses are carried out within a probabilistic framework (see the flowchart for conducting such a RFEM simulation in Figure 3). A linear-elastic, perfectly plastic Tresca soil failure criterion has been used to simulate the undrained clay behaviour. The clay soil parameters are as follows: a unit weight of γs = 20 kN/m3, Young’s modulus of E=100MPa, and Poisson’s ratio of ν=0.3. At the base of the domain, the finite element displacements are fixed, and at the two sides of the domain, it allows only settlements and prevents movements in the perpendicular direction (see Figures 6(a) and 6(b) for the finite element meshes).Figure 6
Finite element meshes and a series of numbered CPT locations for a slope cut problem (a) without and (b) with a foundation layer (dashed lines above the slope face indicate the excavated soil mass, numbers at the bottom of each mesh correspond to Gauss point locations within the corresponding column of finite elements, and an example CPT sample at locationi=0 is indicated).
(a)
(b)Note that it is the2×2 Gauss points within each 8-node finite element that the random field cell values are mapped onto. The random cells can also be mapped onto the finite element levels, but they are not as accurate in simulating the spatial variability [35, 37]. Note also that both the conditional and unconditional random fields of square cells have been mapped onto a finite element mesh with rectangular elements of aspect ratio of 2.0 (Figure 6) to save computational time for the finite element analyses [1, 26, 28]. In this study, the random cell size is 0.125×0.125m, and the random cells are mapped onto the Gauss points within the nonsquare finite elements of size 0.25×0.5m. The readers are referred to [37] for a detailed description of the mapping.


### 4.1. Example 1: A Slope Founded on a Firm Base

The first example considers constructing a 5 m high slope (Figure6(a)), at a proposed slope angle of 45°, in a heterogeneous clay deposit of depth 5 m, overlying on a firm base and whose undrained shear strength cu is spatially varying: it is characterised by a mean of 40 kPa, a standard deviation of 8 kPa, a vertical SOF of θv = 1.0 m, and a horizontal SOF of θh = 6.0, 10, and 14 m for three different cases of horizontal spatial variability, respectively. The coefficient of variation of 0.2 and the chosen values of the scale of fluctuation are well within the recommended ranges for undrained shear strength [35]. This example investigates the effect of the CPT locations on the uncertainty (i.e., in terms of standard deviation) of the realised factor of safety and the effect of CPT distances in the cross section (i.e., y−z plane) during the second phase of sampling.A cross section through the proposed slope is shown in Figure6(a), together with 22 possible locations to site the CPTs (i=0,1,…,21). Note that, in the figure, the finite element size is 0.5 m in the horizontal direction and 0.25 m in the vertical direction. The block of soil is of dimensions 5×11m before the slope is excavated, and the CPTs are carried out before the excavation.The RFEM simulations were carried out both with conditional and unconditional random fields, usingNmc = 300 realisations per Monte Carlo simulation, to investigate the changes in the structure response (i.e., the realised factor of safety) as the sampling (i.e., conditioning) location changes from i=0 to 21. Figures 7(a), 7(c), and 7(e) show that after incorporating the available CPT measurements at i=16, the normalised histograms and fitted probability density distributions of the realised factor of safety for the conditional simulation become narrower, indicating a reduction in the uncertainty of the realised factor of safety after conditioning. Figures 7(b), 7(d), and 7(f) show the example realisations of random fields conditioned upon the CPT at i=16 together with the CPT profiles and a vertical line indicating the CPT location. It is noted that each case of θh corresponds to a different reference “real” field, manifested by the different CPT profiles at i=16. Also, when the CPT profile indicates a strong mean, the updated probability distribution of factor of safety moves to the right of the unconditional distribution. It is also interesting to note that for θh=10m, the PDF moves to the left due to the relatively small values of the CPT profile at the bottom where the failure surface initiates.Figure 7
Simulation results in the first phase of example 1 (based onθv = 1.0 m and 300 realisations per Monte Carlo simulation, CPT location fixed at i=16): (a–b) θh = 6 m; (c–d) θh = 10 m; (e–f) θh = 14 m. (a) PDF of realised factor of safety. (b) Example realisation and CPT. (c) PDF of realised factor of safety. (d) Example realisation and CPT. (e) PDF of realised factor of safety. (f) Example realisation and CPT.
(a)
(b)
(c)
(d)
(e)
(f)Figure8 shows the sampling efficiency indices at various CPT locations, for three different cases with θh=6,10 and 14 m. An optimum position to locate the CPT is clearly seen in the figure; i.e., where the sampling efficiency is a maximum. In this case, i=16, i.e., the optimal CPT location is at the midpoint of the proposed slope. This agrees with the findings in [14]. It is also seen that, when i = 0–1, the improvement is insignificantly small, with Ise close to 1.0. This is due to the various failure surfaces generally passing through zones that are away from the left-hand boundary; at the most right of the mechanistic influence zone, the shear strength is only weakly correlated to values at the sampling location (i.e., i = 0–1). This is especially so for θh=6 and 10 m cases. It is interesting to note that, when i = 21, i.e., the CPT located at the slope toe, the increase in sampling efficiency is still noticeable, that is, sampling at this location is more significant than when sampling at i=0−1. This is due to slope toe area being the most likely zone where slope failure initiates. Therefore, the sampling location relative to the potential failure surfaces matters in the slope stability assessment.Figure 8
Sampling efficiency indices at different sampling locations for various values ofθh (1 CPT sampling), example 1.Note that, for eachθh in the above analysis, the same reference 2D random field is generated first for the rectangular soil block and then used to represent the “real” field variability when generating the conditional random fields in each RFEM analysis. The conditioning step is undertaken first for the rectangular soil block, and the conditional field is then mapped onto the finite element mesh (i.e., the sloped area). This is to ensure that they are consistent with practice, that is, sampling the ground first and then excavating the soil to form a slope. Hence, for CPT locations at i = 12 to 21, although fewer CPT data points are directly used for the elements in the slope face area of the finite element mesh, they nevertheless affect all other cell values (i.e., primarily those cells that are within the correlation distance) through the spatial correlation (i.e., stronger horizontal correlation in the case of anisotropic spatial correlation as is in the current study).Figure8 also shows that, for θh = 10 m and θh = 14 m, the increase in sampling efficiency relative to the unconditional case is greater. This is due to the decreased kriging variance σe2 as θh increases. Moreover, higher values of Ise are obtained for CPT located near the left and right boundaries (especially the right boundary where failure initiates and propagates backwards). Again, this is due to the higher horizontal correlation of soil properties for higher θh values. Also, it is observed that the efficiency of sampling CPT at some distance to the left or right of the midpoint of the slope for higher values of θh can be as much as that of sampling CPT near the slope midpoint for smaller values of θh. For example, the sampling efficiency indices at i = 11 and 20 for θh = 10 m are about the same as that at i = 16 for θh = 6 m.In some cases, a second phase of site investigation may be warranted. This study looks at the influence of a second CPT test (at positionj) at the second stage. The above procedure is repeated by changing j in the range of 0-21 to find the best locations for the new CPT, assuming that the first CPT sample has been taken at its optimal location, i.e., i = 16. Parts of the results are shown in Figure 9 for the cases of θh = 6, 10, and 14 m. Figures 9(a), 9(c), and 9(e) show the normalised histograms and fitted probability density distributions of the realised factor of safety for the conditional simulation of one CPT at i = 16 and for the conditional simulation of an additional CPT at position j = 8. The same horizontal axis scales are used in the three part figures to facilitate comparison with Figures 7(a), 7(c), and 7(e) where the PDF for the unconditional simulation is plotted. It is evident that the PDF has further narrowed down and the confidence in the slope project (either success or failure) has been further improved by the second phase of CPT test. Figures 9(b), 9(d), and 9(f) show the example realisations of random fields conditioned upon 2 CPTs at i=16 and j=8 together with the CPT profiles and two vertical lines indicating the CPT locations. Again, when the 2nd CPT profile indicates a strong mean, the updated probability distribution of factor of safety moves to the right of the distribution based on the 1st CPT only (i.e., for θh = 10 and 14 m). It is noted that, for θh=6m, the PDF moves to the left due to the relatively small mean values of the CPT profile. Should more CPTs be sampled at other places in the cross section, the PDF would gradually narrow down to a definite value, that is, the factor of safety of the reference “real” field.Figure 9
Simulation results in the second phase of example 1 (based onθv = 1.0 m and 300 realisations per Monte Carlo simulation, first phase CPT location fixed at i=16, second phase CPT location j=8): (a–b) θh = 6 m; (c–d) θh = 10 m; (e–f) θh = 14 m. (a) PDF of realised factor of safety. (b) Example realisation and CPT. (c) PDF of realised factor of safety. (d) Example realisation and CPT. (e) PDF of realised factor of safety. (f) Example realisation and CPT.
(a)
(b)
(c)
(d)
(e)
(f)Figure10 shows the sampling efficiency indices for the three cases of θh, with each line representing the influence of various locations j of the second CPT. The figure suggests that the optimal location for conducting the second CPT can be at either side of the slope midpoint (at roughly 4 m distance to the left of the slope midpoint or about 1 m distance to the right of the slope midpoint for the slope geometry considered, i.e., at those locations that equally subdivide the domain to the left or right of the slope midpoint). It is also seen that it is better to take the second CPT at the left side than at the right side. This finding can be largely explained by defining an uncertainty reduction ratio [38] (nyandnz denote the number of target locations or cells to be kriged in the horizontal and vertical directions, respectively),(18)u=∑i=1ny∑j=1nzσei,jnynzσ,and evaluating u for the various sampling strategies of varying second CPT location j. Figure 11 shows the uncertainty ratios as the second CPT column position j changes from 0 to 20 with a step size of 2, and j = 21 for the three cases investigated. It is seen that the locations where the uncertainty ratios reach the minimum roughly coincide with those that are associated with the largest sampling efficiency indices. Also, the uncertainty reduction ratio is smaller on the left side than on the right side. The slight deviation in the optimal locations (to the left) is believed to be attributed to the fact that the uncertainty ratio is defined regardless of the failure influence zone. This highlights the necessity of carrying out a structural response analysis; a kriging error variance analysis alone surely saves the computational expenses associated with the Monte Carlo analysis; however, the mechanical effect is lost when doing so.Figure 10
Sampling efficiency indices at different sampling locations for various values ofθh (2 CPT sampling, 1st CPT fixed at i=16), example 1.Figure 11
Uncertainty ratio for different values ofθh (2 CPT sampling, 1st CPT fixed at i=16), example 1.Conditional simulations considering different numbers of CPTs in the same domain (i.e., a shorter distanceΔ between adjacent profiles with more tests) may be undertaken for different cases of soil spatial variability, to further investigate the effect or impact of CPT intensity on the uncertainty in the slope stability. For example, if budgets allow a third phase site investigation, more CPTs may be sampled halfway in between those at the 2nd phase. The process may be continued until Δ is equal to half the horizontal scale of fluctuation, beyond which little can be gained for the uncertainty reduction [2].


### 4.2. Example 2: A Slope with a Foundation Layer

The second example considers the excavation of a45°, 5 m high slope (Figure 6(b)), in a heterogeneous clay deposit of depth 8 m and of which the undrained shear strength can be statistically characterised by its mean 40 kPa, standard deviation 8 kPa, vertical SOF θv = 1.0 m, and horizontal SOF θh = 10, 14, and 20 m for three different cases of horizontal spatial variability, respectively. Again, the coefficient of variation of 0.2 and the chosen values of the scale of fluctuation are well within the recommended ranges for undrained shear strength [35]. Similar investigations have been carried out in this example to look at the effect of the CPT locations on the uncertainty (i.e., in terms of standard deviation) of the slope stability, and this is then followed by finding the optimal CPT distances in the cross section (i.e., y−z plane).The slope cross section and 40 possible locations to take the CPTs (i=0,1,…,39) is shown in Figure 6(b). Again, the CPTs are carried out before the slope cut excavation, this time in a 8×20m soil deposit (finite element size is 0.5 m in the horizontal direction and 0.25 m in the vertical direction.). The investigation involves both conditional and unconditional RFEM simulations, using Nmc = 300 realisations per simulation. The idea is still to investigate how the slope stability uncertainty (i.e., in terms of the realised factor of safety) is affected by the changes in the sampling location.Similar to Figure8, the sampling efficiency results as a function of the horizontal scale of fluctuation θh and sampling location i is summarised in Figure 12. It is seen that the optimal locations for taking the CPT are at i = 21–24. A higher value of θh results in higher values of the sampling efficiency indices, due to a smaller kriging variance σe2 associated with a higher correlation. It is noted that, for small values of θh, the sampling efficiency indices drop below 1.0 for those locations at the two side boundaries. As mentioned earlier, this is due to the kriging variance being larger than the input shear strength variance (i.e., σe2≥64kPa2) for points far from the measurement points. The optimal location being at or close to the slope crest agrees with [14] for slope stability assessments of slope gradient 3 : 1 (horizontal to vertical), and similar results (though not shown in this paper) have been found for other slope ratios (e.g., 2 : 1) in the current study.Figure 12
Sampling efficiency indices at different sampling locations for various values ofθh (1 CPT sampling), example 2.To further verify the above finding regarding the best location to sample the 2nd CPT in the 2nd phase of site investigation, 39 sets of conditional simulations for differing 2nd phase sampling locationsj while the first CPT was fixed at i=21 have been carried out. The results are shown in Figure 13 for the three values of θh; it is seen that it is indeed the best to sample halfway between the 1st CPT and the side boundaries, within the influence zone of the failure surface. The same reason as in example 1 applies here; that is, the uncertainty ratio (equation (18)) reaches a minimum at those locations.Figure 13
Sampling efficiency indices at different sampling locations for various values ofθh (2 CPT sampling, 1st CPT fixed at i=21), example 2.


## 4.1. Example 1: A Slope Founded on a Firm Base

The first example considers constructing a 5 m high slope (Figure6(a)), at a proposed slope angle of 45°, in a heterogeneous clay deposit of depth 5 m, overlying on a firm base and whose undrained shear strength cu is spatially varying: it is characterised by a mean of 40 kPa, a standard deviation of 8 kPa, a vertical SOF of θv = 1.0 m, and a horizontal SOF of θh = 6.0, 10, and 14 m for three different cases of horizontal spatial variability, respectively. The coefficient of variation of 0.2 and the chosen values of the scale of fluctuation are well within the recommended ranges for undrained shear strength [35]. This example investigates the effect of the CPT locations on the uncertainty (i.e., in terms of standard deviation) of the realised factor of safety and the effect of CPT distances in the cross section (i.e., y−z plane) during the second phase of sampling.A cross section through the proposed slope is shown in Figure6(a), together with 22 possible locations to site the CPTs (i=0,1,…,21). Note that, in the figure, the finite element size is 0.5 m in the horizontal direction and 0.25 m in the vertical direction. The block of soil is of dimensions 5×11m before the slope is excavated, and the CPTs are carried out before the excavation.The RFEM simulations were carried out both with conditional and unconditional random fields, usingNmc = 300 realisations per Monte Carlo simulation, to investigate the changes in the structure response (i.e., the realised factor of safety) as the sampling (i.e., conditioning) location changes from i=0 to 21. Figures 7(a), 7(c), and 7(e) show that after incorporating the available CPT measurements at i=16, the normalised histograms and fitted probability density distributions of the realised factor of safety for the conditional simulation become narrower, indicating a reduction in the uncertainty of the realised factor of safety after conditioning. Figures 7(b), 7(d), and 7(f) show the example realisations of random fields conditioned upon the CPT at i=16 together with the CPT profiles and a vertical line indicating the CPT location. It is noted that each case of θh corresponds to a different reference “real” field, manifested by the different CPT profiles at i=16. Also, when the CPT profile indicates a strong mean, the updated probability distribution of factor of safety moves to the right of the unconditional distribution. It is also interesting to note that for θh=10m, the PDF moves to the left due to the relatively small values of the CPT profile at the bottom where the failure surface initiates.Figure 7
Simulation results in the first phase of example 1 (based onθv = 1.0 m and 300 realisations per Monte Carlo simulation, CPT location fixed at i=16): (a–b) θh = 6 m; (c–d) θh = 10 m; (e–f) θh = 14 m. (a) PDF of realised factor of safety. (b) Example realisation and CPT. (c) PDF of realised factor of safety. (d) Example realisation and CPT. (e) PDF of realised factor of safety. (f) Example realisation and CPT.
(a)
(b)
(c)
(d)
(e)
(f)Figure8 shows the sampling efficiency indices at various CPT locations, for three different cases with θh=6,10 and 14 m. An optimum position to locate the CPT is clearly seen in the figure; i.e., where the sampling efficiency is a maximum. In this case, i=16, i.e., the optimal CPT location is at the midpoint of the proposed slope. This agrees with the findings in [14]. It is also seen that, when i = 0–1, the improvement is insignificantly small, with Ise close to 1.0. This is due to the various failure surfaces generally passing through zones that are away from the left-hand boundary; at the most right of the mechanistic influence zone, the shear strength is only weakly correlated to values at the sampling location (i.e., i = 0–1). This is especially so for θh=6 and 10 m cases. It is interesting to note that, when i = 21, i.e., the CPT located at the slope toe, the increase in sampling efficiency is still noticeable, that is, sampling at this location is more significant than when sampling at i=0−1. This is due to slope toe area being the most likely zone where slope failure initiates. Therefore, the sampling location relative to the potential failure surfaces matters in the slope stability assessment.Figure 8
Sampling efficiency indices at different sampling locations for various values ofθh (1 CPT sampling), example 1.Note that, for eachθh in the above analysis, the same reference 2D random field is generated first for the rectangular soil block and then used to represent the “real” field variability when generating the conditional random fields in each RFEM analysis. The conditioning step is undertaken first for the rectangular soil block, and the conditional field is then mapped onto the finite element mesh (i.e., the sloped area). This is to ensure that they are consistent with practice, that is, sampling the ground first and then excavating the soil to form a slope. Hence, for CPT locations at i = 12 to 21, although fewer CPT data points are directly used for the elements in the slope face area of the finite element mesh, they nevertheless affect all other cell values (i.e., primarily those cells that are within the correlation distance) through the spatial correlation (i.e., stronger horizontal correlation in the case of anisotropic spatial correlation as is in the current study).Figure8 also shows that, for θh = 10 m and θh = 14 m, the increase in sampling efficiency relative to the unconditional case is greater. This is due to the decreased kriging variance σe2 as θh increases. Moreover, higher values of Ise are obtained for CPT located near the left and right boundaries (especially the right boundary where failure initiates and propagates backwards). Again, this is due to the higher horizontal correlation of soil properties for higher θh values. Also, it is observed that the efficiency of sampling CPT at some distance to the left or right of the midpoint of the slope for higher values of θh can be as much as that of sampling CPT near the slope midpoint for smaller values of θh. For example, the sampling efficiency indices at i = 11 and 20 for θh = 10 m are about the same as that at i = 16 for θh = 6 m.In some cases, a second phase of site investigation may be warranted. This study looks at the influence of a second CPT test (at positionj) at the second stage. The above procedure is repeated by changing j in the range of 0-21 to find the best locations for the new CPT, assuming that the first CPT sample has been taken at its optimal location, i.e., i = 16. Parts of the results are shown in Figure 9 for the cases of θh = 6, 10, and 14 m. Figures 9(a), 9(c), and 9(e) show the normalised histograms and fitted probability density distributions of the realised factor of safety for the conditional simulation of one CPT at i = 16 and for the conditional simulation of an additional CPT at position j = 8. The same horizontal axis scales are used in the three part figures to facilitate comparison with Figures 7(a), 7(c), and 7(e) where the PDF for the unconditional simulation is plotted. It is evident that the PDF has further narrowed down and the confidence in the slope project (either success or failure) has been further improved by the second phase of CPT test. Figures 9(b), 9(d), and 9(f) show the example realisations of random fields conditioned upon 2 CPTs at i=16 and j=8 together with the CPT profiles and two vertical lines indicating the CPT locations. Again, when the 2nd CPT profile indicates a strong mean, the updated probability distribution of factor of safety moves to the right of the distribution based on the 1st CPT only (i.e., for θh = 10 and 14 m). It is noted that, for θh=6m, the PDF moves to the left due to the relatively small mean values of the CPT profile. Should more CPTs be sampled at other places in the cross section, the PDF would gradually narrow down to a definite value, that is, the factor of safety of the reference “real” field.Figure 9
Simulation results in the second phase of example 1 (based onθv = 1.0 m and 300 realisations per Monte Carlo simulation, first phase CPT location fixed at i=16, second phase CPT location j=8): (a–b) θh = 6 m; (c–d) θh = 10 m; (e–f) θh = 14 m. (a) PDF of realised factor of safety. (b) Example realisation and CPT. (c) PDF of realised factor of safety. (d) Example realisation and CPT. (e) PDF of realised factor of safety. (f) Example realisation and CPT.
(a)
(b)
(c)
(d)
(e)
(f)Figure10 shows the sampling efficiency indices for the three cases of θh, with each line representing the influence of various locations j of the second CPT. The figure suggests that the optimal location for conducting the second CPT can be at either side of the slope midpoint (at roughly 4 m distance to the left of the slope midpoint or about 1 m distance to the right of the slope midpoint for the slope geometry considered, i.e., at those locations that equally subdivide the domain to the left or right of the slope midpoint). It is also seen that it is better to take the second CPT at the left side than at the right side. This finding can be largely explained by defining an uncertainty reduction ratio [38] (nyandnz denote the number of target locations or cells to be kriged in the horizontal and vertical directions, respectively),(18)u=∑i=1ny∑j=1nzσei,jnynzσ,and evaluating u for the various sampling strategies of varying second CPT location j. Figure 11 shows the uncertainty ratios as the second CPT column position j changes from 0 to 20 with a step size of 2, and j = 21 for the three cases investigated. It is seen that the locations where the uncertainty ratios reach the minimum roughly coincide with those that are associated with the largest sampling efficiency indices. Also, the uncertainty reduction ratio is smaller on the left side than on the right side. The slight deviation in the optimal locations (to the left) is believed to be attributed to the fact that the uncertainty ratio is defined regardless of the failure influence zone. This highlights the necessity of carrying out a structural response analysis; a kriging error variance analysis alone surely saves the computational expenses associated with the Monte Carlo analysis; however, the mechanical effect is lost when doing so.Figure 10
Sampling efficiency indices at different sampling locations for various values ofθh (2 CPT sampling, 1st CPT fixed at i=16), example 1.Figure 11
Uncertainty ratio for different values ofθh (2 CPT sampling, 1st CPT fixed at i=16), example 1.Conditional simulations considering different numbers of CPTs in the same domain (i.e., a shorter distanceΔ between adjacent profiles with more tests) may be undertaken for different cases of soil spatial variability, to further investigate the effect or impact of CPT intensity on the uncertainty in the slope stability. For example, if budgets allow a third phase site investigation, more CPTs may be sampled halfway in between those at the 2nd phase. The process may be continued until Δ is equal to half the horizontal scale of fluctuation, beyond which little can be gained for the uncertainty reduction [2].


## 4.2. Example 2: A Slope with a Foundation Layer

The second example considers the excavation of a45°, 5 m high slope (Figure 6(b)), in a heterogeneous clay deposit of depth 8 m and of which the undrained shear strength can be statistically characterised by its mean 40 kPa, standard deviation 8 kPa, vertical SOF θv = 1.0 m, and horizontal SOF θh = 10, 14, and 20 m for three different cases of horizontal spatial variability, respectively. Again, the coefficient of variation of 0.2 and the chosen values of the scale of fluctuation are well within the recommended ranges for undrained shear strength [35]. Similar investigations have been carried out in this example to look at the effect of the CPT locations on the uncertainty (i.e., in terms of standard deviation) of the slope stability, and this is then followed by finding the optimal CPT distances in the cross section (i.e., y−z plane).The slope cross section and 40 possible locations to take the CPTs (i=0,1,…,39) is shown in Figure 6(b). Again, the CPTs are carried out before the slope cut excavation, this time in a 8×20m soil deposit (finite element size is 0.5 m in the horizontal direction and 0.25 m in the vertical direction.). The investigation involves both conditional and unconditional RFEM simulations, using Nmc = 300 realisations per simulation. The idea is still to investigate how the slope stability uncertainty (i.e., in terms of the realised factor of safety) is affected by the changes in the sampling location.Similar to Figure8, the sampling efficiency results as a function of the horizontal scale of fluctuation θh and sampling location i is summarised in Figure 12. It is seen that the optimal locations for taking the CPT are at i = 21–24. A higher value of θh results in higher values of the sampling efficiency indices, due to a smaller kriging variance σe2 associated with a higher correlation. It is noted that, for small values of θh, the sampling efficiency indices drop below 1.0 for those locations at the two side boundaries. As mentioned earlier, this is due to the kriging variance being larger than the input shear strength variance (i.e., σe2≥64kPa2) for points far from the measurement points. The optimal location being at or close to the slope crest agrees with [14] for slope stability assessments of slope gradient 3 : 1 (horizontal to vertical), and similar results (though not shown in this paper) have been found for other slope ratios (e.g., 2 : 1) in the current study.Figure 12
Sampling efficiency indices at different sampling locations for various values ofθh (1 CPT sampling), example 2.To further verify the above finding regarding the best location to sample the 2nd CPT in the 2nd phase of site investigation, 39 sets of conditional simulations for differing 2nd phase sampling locationsj while the first CPT was fixed at i=21 have been carried out. The results are shown in Figure 13 for the three values of θh; it is seen that it is indeed the best to sample halfway between the 1st CPT and the side boundaries, within the influence zone of the failure surface. The same reason as in example 1 applies here; that is, the uncertainty ratio (equation (18)) reaches a minimum at those locations.Figure 13
Sampling efficiency indices at different sampling locations for various values ofθh (2 CPT sampling, 1st CPT fixed at i=21), example 2.


## 5. Conclusions

Conditional simulation can be used as a useful aid to geotechnical sampling design. Random fields conditioned upon CPT measurements are used in this paper to investigate the influence of CPT location and distance on slope stability uncertainty. The potential use of conditional simulation in geotechnical site exploration and cost-effective designs is illustrated through two numerical slope cut examples. The results showed that the confidence in the stability of a slope can be increased by conditional simulations honoring CPT measurement data. Indeed, the unconditional slope stability simulation based on the statistics of soil properties only (without explicit consideration of the constraints of CPT measurements at the sampling locations) results in a “prior” distribution of the factor of safety. In contrast, after the inclusion of the “real” spatial distribution of all CPTs, a “posterior” distribution of the structure performance (i.e., factor of safety in this paper) can be obtained. Therefore, the probability of failure (or reliability) from a conditional slope stability simulation can be viewed as a conditional probability of failure (or reliability). The updating and improvement of the probability density distributions of factor of safety in the two numerical examples clearly demonstrate this. For both slope examples, with or without a foundation layer, the conditional slope stability simulation led to a more cost-effective stability design (i.e., “effective” in the sense that the uncertainty in the response reduces) through eliminating unrealistic soil property realisations that are not consistent with measured CPT profiles.The study also shows the updating of the response probability density function in a 2nd phase of site investigation and the improvement of the confidence in the probability of failure or survival in the 2nd phase. In fact, a site investigation may be carried out in multiple stages in many cases, with the initial analysis results (e.g., conditional simulation based on the 1st CPT) guiding further field tests (e.g., the 2nd CPT location or its distance from the 1st one). As demonstrated in the two examples, in the case of two-phase site investigation, it is possible to find out the optimum distance between the 2nd CPT and the 1st CPT, i.e., the optimal location for the additional testing in the 2nd phase. The potential use of the method in directing site exploration programs is highlighted and the efficient utilisation of field measurements is emphasized.For the examples considered in this paper, the optimal sampling location lies near the slope crest to the midpoint of the slope, and the optimal sampling distance between CPTs for the 2nd phase was identified to equally subdivide the zones to the left or right of the 1st test. The general suggestion based on the current investigation is to sample at the best location in the first phase and then sample to either side of the first test at some locations that equally subdivide the mechanical influence domain to the left or right of the first test. This subdividing process goes on until the distance between CPTs is equal to or just less than half the horizontal scale of fluctuation [2]. The implication is that some prior knowledge should be obtained first about the mechanical influence zone for the problem at hand, perhaps in a deterministic context. It should again be mentioned that the current study is based on a 2D cross section. In realistic 3D applications with 3D spatial variability, the first phase of investigation often involves a number of CPTs in the third dimension. This first row of CPTs (i.e., in the 3rd dimension) serves as a first data set to be used to estimate the horizontal scale of fluctuation, and this first estimate should guide the second row of CPTs in the cross section (i.e., relative to the first row) in the 2nd phase of site investigation (i.e., based on the findings in this paper and the previous findings in 3D and a reasonable assumption of an isotropic correlation structure in the horizontal depositional plane).



---

*Source: 8267601-2019-09-10.xml*