Delft University of Technology
Uncertainty reduction and sampling efficiency in slope designs using 3D conditional
random fields
Li, Y.; Hicks, M. A.; Vardon, P. J. DOI
10.1016/j.compgeo.2016.05.027 Publication date
2016
Document Version
Accepted author manuscript Published in
Computers and Geotechnics
Citation (APA)
Li, Y., Hicks, M. A., & Vardon, P. J. (2016). Uncertainty reduction and sampling efficiency in slope designs using 3D conditional random fields. Computers and Geotechnics, 79, 159-172.
https://doi.org/10.1016/j.compgeo.2016.05.027 Important note
To cite this publication, please use the final published version (if applicable). Please check the document version above.
Copyright
Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy
Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.
This work is downloaded from Delft University of Technology.
Uncertainty reduction and sampling efficiency in slope designs using 3D
1
conditional random fields
2
Y.J. Li, M.A. Hicks, P.J. Vardon
3
Section of Geo-Engineering, Department of Geoscience & Engineering, Faculty of Civil Engineering
4
& Geosciences, Delft University of Technology, Delft, The Netherlands
5
Email: y.li-4@tudelft.nl; m.a.hicks@tudelft.nl; p.j.vardon@tudelft.nl
6
Cite as:
7
Li, Y. J., Hicks, M. A., & Vardon, P. J. (2016). Uncertainty reduction and sampling
8
efficiency in slope designs using 3D conditional random fields. Computers and
9
Geotechnics, 79: 159-172.
doi:10.1016/j.compgeo.2016.05.02710
free link until 16 August, 2016. http://authors.elsevier.com/a/1THNi,63b~copS
11This is the Authors’ final manuscript. (page number 1-27) 12
Abstract
13A method of combining 3D Kriging for geotechnical sampling schemes with an existing random field 14
generator is presented and validated. Conditional random fields of soil heterogeneity are then linked 15
with finite elements, within a Monte Carlo framework, to investigate optimum sampling locations and 16
the cost-effective design of a slope. The results clearly demonstrate the potential of 3D conditional 17
simulation in directing exploration programmes and designing cost saving structures; that is, by 18
reducing uncertainty and improving the confidence in a project’s success. Moreover, for the problems 19
analysed, an optimal sampling distance of half the horizontal scale of fluctuation was identified. 20
Key words: conditional random fields, Kriging, reliability, sampling efficiency, spatial variability, 21
uncertainty reduction 22
1. Introduction
23Soil properties exhibit three dimensional spatial variability (i.e. heterogeneity). In geotechnical 24
engineering, a site investigation may be carried out, and the data collected and processed in a 25
statistical way to characterise the variability [1–10]. The outcomes of the statistical treatment, e.g. the 26
mean property value, the standard deviation or coefficient of variation, and the spatial correlation 27
distance, may be used as input to a geotechnical model capable of dealing with the spatial variation 28
(e.g. a random field simulation). However, when it comes to making use of the field data, there arises 29
the question: How can we make best use of the available data? The idea is to use the data more 30
effectively, so that it is worth the effort or cost spent in carrying out the investigation, as well as the 31
additional effort in post-processing the data. The aim of this paper is to contribute towards answering 32
this question. 33
For example, cone penetration tests (CPTs) are often carried out in geotechnical field investigations, in 34
order to obtain data used in implementing the design of a structure. The amount of data from CPT 35
measurements is often larger than from conventional laboratory tests. This is useful, as a large 36
database is needed to accurately estimate the spatial correlation structure of a soil property. For 37
example, Fenton [3] used a database of CPT profiles from Oslo to estimate the correlation statistics in 38
the vertical direction, and Jaksa et al. [5] used a database from Adelaide to estimate the correlation 39
distances in both the vertical and horizontal directions. 40
In geotechnical engineering, a substantial amount of numerical work has been done using idealised 2D 41
simulations based on collected in-situ data (e.g. [4]), although a 3D simulation would be preferable 42
due to site data generally being collected from a 3D space. However, there are relatively few studies 43
simulating the effect of 3D heterogeneity due to the high computational requirements. Examples 44
include the effect of heterogeneity on shallow foundation settlement [11–13], on steady state seepage 45
[14–16], on seismic liquefaction [17] and on slope reliability [18–27]. 46
The above investigations all used random fields to represent the soil spatial variability and the finite 47
element method to analyse geotechnical performance within a Monte Carlo framework, a form of 48
analysis sometimes referred to as the random finite element method (RFEM) [28]. However, they did 49
not make use of the spatial distribution of related measurement data to constrain the random fields. In 50
other words, for those applications that are based on real field data, many realisations not complying 51
with the field data at the measurement locations will be included in the simulation, which, in turn, will 52
result in an exaggerated range of responses in the analysis of geotechnical performance. 53
Studies on conditional simulations are available in geostatistics in the field of reservoir engineering 54
[29]. However, there are not many studies dealing with soil spatial variability in geotechnical 55
engineering that utilise conditional simulation (some 2D exceptions include, e.g., [6, 30–32]). This is 56
partly due to the smaller amount of data generally available in geotechnical engineering, and partly 57
due to there often not being a computer program specially implemented for those situations where 58
there are sufficient data (e.g. CPT, vane shear test (VST)), especially in 3D. However, unconditional 59
random fields can easily be conditioned to the known measurements by Kriging [29, 33]. Hence, 60
following the previous 2D work of Van den Eijnden and Hicks [31] and Lloret-Cabot et al. [30], this 61
paper seeks to implement and apply conditional simulation in three dimensional space, in order to 62
reduce uncertainty in the field where CPT measurements are carried out. 63
Usually, site investigation plans are designed to follow some regular pattern. For example, a 64
systematic grid of sample locations is generally used, due to its simplicity to implement [5]. Moreover, 65
although there are various sampling plans in terms of layouts, it is found that systematically ordered 66
spatial samples are superior in terms of the quality of estimates at unsampled locations [34]. Therefore, 67
this paper will be devoted to implementing a 3D Kriging algorithm for sampling schemes following a 68
regular grid. This will then be combined with an existing 3D random field generator to implement a 69
conditional simulator. However, extension to irregular sampling patterns is straightforward based on 70
the presented framework. 71
The implemented approach has been applied to two idealised slope stability examples. The first 72
demonstrates how the approach may be used to identify the best locations to conduct borehole testing, 73
and thereby allow an increased confidence in a project’s success or failure to be obtained. While it is 74
very important to pay sufficient attention to the required intensity of a site investigation (i.e. the 75
optimal number of boreholes) with respect to the site-specific spatial variability, as highlighted by 76
Jaksa et al. [12], the first example starts by focusing on the optimum locations for carrying out site 77
investigations for a given number of boreholes, before moving on to consider the intensity of testing. 78
The second example compares different candidate slope designs, in order to choose the best (most 79
cost-effective) design satisfying the reliability requirements. 80
For simplicity, this paper focuses on applications involving only a single soil layer (i.e. a single layer 81
characterised by a statistically homogeneous undrained shear strength), although the extension to 82
multiple soil layers is straightforward. Moreover, the effect of random variation in the boundary 83
locations between different soil layers can also be easily incorporated by conditioning to known 84
boundary locations (e.g. corresponding to where the CPTs have been carried out). 85
2. Theory and Implementation
862.1 Conditioning
87
A conditional random field, which preserves the known values at the measurement locations, can be 88
formed from three different fields [28, 35–36]: 89
(
)
( )
( )
( )
( )
rc ru km ksZ
x
=
Z
x
+
Z
x
−
Z
x
(1) 90wherexdenotes a location in space,Zrc( )x is the conditionally simulated random field,Zru( )x is the 91
unconditional random field,Zkm( )x is the Kriged field based on measured values at xi(i=1, 2,2,N), 92
( ) ks
Z x is the Kriged field based on unconditionally (or randomly) simulated values at the same 93
positions xi(i=1, 2,2,N), and N is the number of measurement locations. 94
The unconditional random field can be simulated via several methods [37]; for example, interpolated 95
autocorrelation [38], covariance matrix decomposition, discrete Fourier transform or Fast Fourier 96
transform, turning bands, local average subdivision (LAS), and Karhunen–Loeve expansion [39], 97
among others. The LAS method [40] is used in this paper. The Kriged fields are obtained by Kriging 98
[41], which has found extensive usage in geostatistics [42–43]. The LAS and Kriging methods are 99
briefly reviewed in the following sections. 100
2.2 Anisotropic random field generation using 3D LAS
101
The LAS method [40, 44] is used herein to generate the unconditional random fields, using statistics 102
(i.e. mean, variance and correlation structure) based on the observed field data. The LAS method 103
proceeds in a recursive fashion, by progressively subdividing the initial domain into smaller cells, until 104
the random process is represented by a series of local averages. The major advantage is its ability to 105
produce random fields of local averages whose statistics are consistent with the field resolution; that is, 106
it maintains a constant mean over all levels of subdivision, and ensures reduced variances as a function 107
of cell size based on variance reduction theory [45], taking account of spatial correlations between 108
local averages within each level and across levels. 109
The following covariance function is used in the subdivision process: 110
( )
(
)
2 1 2 2 3 2 1 2 3 1 2 32
2
2
,
,
exp
C
C
τ τ τ
σ
τ
τ
τ
θ
θ
θ
=
=
−
−
+
τ
(2) 111where
σ
2is the variance of the soil property,τ
is the lag vector, and θ1, θ2 and θ3, and τ1, τ2 and 112τ3 are the respective scales of fluctuation and lag distances in the vertical and two lateral coordinate 113
directions, respectively. Herein, an isotropic random field is initially generated by setting θ1 = θ2 = θ3 114
= θiso; i.e. so that θiso equals the horizontal scale of fluctuation, θh. This field is then squashed in the 115
vertical direction to give the target vertical scale of fluctuation, θv. The 3D LAS implementation of 116
Spencer [25] has been used in this paper, and the reader is referred to Spencer [25] and Hicks and 117
Spencer [19] for more details. Note also that a truncated normal distribution has been used to describe 118
the pointwise variation in material properties [19]. 119
2.3 Kriging
120
In contrast to conventional deterministic interpolation techniques, such as moving least squares and 121
the radial point interpolation method, Kriging incorporates the variogram (or covariance) into the 122
interpolation procedure; specifically, information on the spatial correlation of the measured points is 123
used to calculate the weights. Moreover, standard errors of the estimation can also be obtained, 124
indicating the reliability of the estimation and the accuracy of the prediction. Kriging is a method of 125
interpolation for which the interpolated values are modelled by a Gaussian process governed by prior 126
covariances and for which confidence intervals can be derived. While interpolation methods based on 127
other criteria need not yield the most likely intermediate values, Kriging provides a best linear 128
unbiased prediction of the soil properties (Z) between known data [43, 46] by assuming the stationarity 129
of the mean and of the spatial covariances, or variograms. A brief review is first given to facilitate 130
understanding of the implementation. 131
Suppose that Z Z1, 2,2,ZNare observations of the random field
Z x
( )
at points x x1, 2,2,xN(i.e. 132( ) ( 1, 2, , )
i i
Z =Z x i= 2 N ). The best linear unbiased estimation (i.e. Zˆ) of the soil property at some 133 locationx0is given by 134 0 0 1 1
ˆ ( )
(
) ( )
N N i i i i i iZ
λ
Z
λ
Z
= ==
∑
=
∑
x
x
x
(3) 135in which N denotes the total number of observations and
λ
idenotes the unknown weighting factor 136associated with observation point xi, which needs to be determined. 137
The weights in equation (3), for the estimation at any location
x
0, can be found by minimising the 138variance (
σ
e2) of the Kriging error Zˆ−Z , which is given as 139(
)
(
)
(
)
(
)
(
(
)
)
(
(
)
)
(
)
(
) (
)
2 2 2 2 2 0 00 1 1 1 2 2 2 0 0 0 1 1 1 0 0 0 1 1 1ˆ
ˆ
var
(
) 2
(
) (
)
2
2
N N N e i j ij i i i j i N N N i j i j i i i j i N N N i j i j i i i j iZ
Z
E
Z
Z
c
c
c
C
C
C
σ
λ λ
σ
λ
σ
σ
λ λ
σ
λ
σ
σ
λ λ γ
λ γ
γ
= = = = = = = = =
=
−
=
−
=
−
−
−
+
−
=
−
−
−
−
−
+
−
−
= −
−
+
−
−
−
∑∑
∑
∑∑
∑
∑∑
∑
x
x
x
x
x
x
x
x
x
x
x
x
(4) 140 6where var() denotes the variance operator and E[] is the expectation operator, cij =C
(
xi−xj)
is the 141covariance between ( )Z x and i Z x( j),ci0=C
(
xi−x0)
is the covariance between ( )Z x and i Z x( 0), 142and c00=C
(
x0−x0)
=C(0) =σ2 is the variance of 0( )
Z x , which is estimated at the target location 143
0
x . The rearrangement in equation (4) makes use of the relationship between a variogram
γ
( )
τ and a 144covariance function C
( )
τ ( i.e.γ
( )
τ =C( )
0 −C( )
τ =σ
2−C( )
τ ) and the condition1 1 N i i λ = =
∑
(in 145order to ensure that the estimator is unbiased, i.e. E Z
(
ˆ−Z)
= , the weights must sum to one). 0 146To minimise the error variance (i.e. equation (4)), the Lagrange method is used [43]. The weights can 147
then be found by solving the system of equations of size N + 1, for a constant mean: 148 1 1 1 1 1 0 1 0
(
)
(
)
1
(
)
(
)
(
)
1
(
)
1
1
0
1
N N N N N Nγ
γ
λ
γ
γ
γ
λ
γ
µ
−
−
−
=
−
−
−
x
x
x
x
x
x
x
x
x
x
x
x
(5) 149in which
µ
is the Lagrangian parameter. For a mean following some trend, the modification to 150equation (5) is straightforward and interested readers are referred to Fenton [46]. 151
Equation (5) may be expressed as 152
lhs x = rhs
γ λ γ (6) 153
Once equation (5) is solved, the estimated error variance can be expressed by 154
(
)
(
)
T 2 1 0 1 N e i i lhs rhs rhs iσ
µ
λ γ
− == +
∑
x
−
x
=
γ
γ
γ
(7) 155where
λ
iis a function of the relative positioning of pointsxiand x0. 156Note that the left-hand-side matrix γlhsis a function of only the observation point locations and 157
covariance between them. Therefore, it only needs to be inverted once, and then equations (5) and (3) 158
used repeatedly for building up the field of best estimates at different locations in space. In contrast, 159
the right-hand-side vector γrhschanges as a function of the spatial point x0, resulting in different 160
weight vectors λxthat are used in equation (3) to get the estimates (point by point) in the domain of 161
interest. 162
In geotechnical engineering, a sampling strategy following some pattern is generally adopted [1]. For 163
example, CPT sampling is often planned in the form of a regular grid on the ground surface [5]. It is 164
therefore desirable to implement the above Kriging algorithm in the context of some sampling design 165
with a regular pattern. While it is straightforward to implement in 2D, it is less so when implemented 166
in 3D. The most fundamental part is how the left-hand-side matrix of equation (6) is formed. The 167
authors have implemented 3D Kriging for the regular grid sampling strategy shown in Figure 1. The 168
way to set up the left-hand-side matrix and the right-hand-side vector is presented in the Appendix. 169
2.4 Computational efficiency
170
There are two aspects involved in the computational efficiency of the above Kriging implementation. 171
One is the total number of equations, which depends on the total number of data points (N = k×m×n, 172
where k and m are the number of CPT rows in the x and y directions respectively, and n is the number 173
of data points for each CPT profile, see Figure 1) contributing to the left-hand-side matrix; the other is 174
the number of points in the field (nf = nx×ny×nz, where nx, ny and nz are the number of points in the 175
three Cartesian directions) that need to be Kriged (i.e. how many times the algorithm will need to be 176
repeated, except for inverting the left-hand-side matrix). The higher the required field resolution (nf) 177
and the greater the total number of known data points (N), the longer the Kriging will take. In the case 178
of the CPT arrangement in Fig. 1, the size of matrix γlhs (see equation (5) or equation (A1)) and the 179
size of vector γrhs(see equation (5) or equation (A3)) in 3D are m2 and m times larger than those in 2D 180
(i.e. a cross-section in the x–z plane) respectively. The time it takes to Krige a full 3D field depends on 181
the processing time of each individual step and the number of times each step has to be performed. To 182
Krige a field of size nf, conditional to N measurement points, the total time may be approximated by 183
2 3
1 2
( , f) f
t N n ≈c n N +c N (8) 184
where the first term represents the time needed for solving the system of equations for all field points 185
(i.e. nf times) (O(N 2
)) and the second represents the time needed for inverting the matrix
γ
lhs (i.e. only 186once) (O(N3)). The constants c1 and c2 are functions of the CPU speed and the operation, and in this 187
case are in a ratio of approximately 4:1. Additionally, in all practical cases, nf >> N, so that the 188
calculation time depends mainly on the first term in the above equation; that is, on the number of times 189
(i.e. nf times) that the matrix–vector multiplication operation,
1 x lhs rhs −
=
λ
γ
γ
, needs to be performed. 190Note that nf = nx×ny×nz in 3D is ny times nf = nx×nz in 2D and N = k×m×n in 3D is m times N = k×n in 191
2D. In the examples reported in Section 4, all problems investigated are very long in the third 192
dimension compared to the cross-section. Therefore, the time consumed in a 3D analysis is 193
theoretically ny×m 2
times that of a 2D analysis, when neglecting the relatively fast, one-off matrix 194
inversion operation and other computation overheads, such as reading/writing and memory operations. 195
However, despite the significantly greater run-time requirements for Kriging in 3D (as compared to 196
2D), it is still far less than the time consumed in a nonlinear finite element analysis where plasticity 197
iterations are needed. For Example 1 in Section 4, where nx = 20, ny = 100 and nz = 20, it took, in serial 198
and on average, 134 hours in total for the 500 finite element analyses forming each Monte Carlo 199
simulation (3.0 GHz CPU), whereas Kriging 500 times took about 2.4 hours. In contrast, 500 Kriging 200
interpolations for a 2D cross-section analysis took approximately 8.5 seconds. It is noted that the 201
computation time used for Kriging 500 times is significant in comparison with a single finite element 202
analysis, and therefore should not be considered a pre-processing step if utilising parallel computation 203
for the finite element analyses. Therefore, the computing strategy developed to carry out the analyses 204
for Examples 1 and 2 in Section 4 (comprising around 30,000 realisations in total, and involving 205
30,000 3D Kriging interpolations) was to run the analyses in parallel (each Kriging and finite element 206
analysis serially on a single computation node) on the Dutch national grid e-infrastructure with high 207
performance computing clusters. 208
Note that it is possible to prescribe an appropriate neighbourhood size in the algorithm to reduce the 209
computational burden for 3D Kriging. For example, a neighbourhood size of 5×7×n may be used to 210
construct the left-hand-side matrix (see the neighbourhood denoted as a rectangle in Figure A.1(a), i.e. 211
by using only the nearest 4 CPT profiles). That is, only those CPT profiles that have a significant 212
influence (i.e. a lag distance within the range of the scale of fluctuation in equation (2)) on the point to 213
be estimated are used to construct the left-hand-side (LHS) matrix. However, using this strategy, for 214
each point (or each subset of points) to be estimated, the left-hand-side matrix is different and will 215
need to be inverted accordingly, so this could increase the computational time if there are a large 216
number of points or cells to be estimated. Therefore, a choice has to be made, to make sure that the 217
time saved by inverting a smaller matrix, instead of a bigger one, outweighs the time consumed by 218
inverting the left-hand-side matrices for all the (subgroups of) cells to be estimated for the case in 219
which a neighbourhood is used. And, of course, there is a trade-off between the estimation accuracy 220
and time saved when such a neighbourhood approach is used. The accuracy will increase as more 221
available data are used to do the Kriging estimation, and so the neighbourhood size depends on the 222
required accuracy and the scales of fluctuation. 223
Due to the relatively fast inversion of the LHS matrix in the current investigation (the maximum size 224
investigated is N = 500), all CPT profiles have been used for the Kriging in the examples in Section 4. 225
However, one neighbourhood strategy was investigated by using the 4 nearest CPT profiles, and the 226
following uncertainty reduction ratio (a 3D extension to the 1D definition in [31]) has been used to 227
assess the approximation error: 228
(
)
1 1 1 , , y x n z n n e i j k x y z i j k u n n n σ σ = = = =∑∑∑
(9) 229The approximation error may be evaluated by 230 n a u a u u E u − = (10) 231 10
where un and ua are the uncertainty reduction ratios when using a neighbourhood and when all CPT 232
profiles have been used, respectively. 233
One of the sampling strategies from Example 1 (Section 4, Fig. 9(b)) was used to evaluate the 234
approximation error and the results are listed in Table 1. It can be seen that using a neighbourhood of 235
the 4 nearest CPT profiles has been sufficient in this case. 236
3. Validation
237The conditional simulation of a 5 m high (z), 5 m wide (x) and 25 m long (y) clay block, characterised 238
by a spatially varying undrained shear strength, is presented in this section to demonstrate the 239
procedure and the validity of the implementation described in Sections 2.1–2.3 (and the Appendix). 240
The idea is to show how the measured values are honoured, and to check whether or not the statistical 241
properties (e.g. covariance) of the random fields are maintained after conditioning. 242
The block is discretised into 20×100×20 cubic cells, with each cell of dimension 0.25 m. The mean of 243
the undrained shear strength is 40 kPa, and the standard deviation is 8 kPa. The degree of anisotropy 244
of the heterogeneity is ξ = 3, in which ξ = θh/θv and θv = 1.0 m. Five CPT measurement locations in 245
the y direction (at x = 2.5 m) are available, each comprising n = 20 data points at 0.25 m spacing in the 246
vertical direction. These ‘measured’ data have been obtained by sampling from a single independent 247
realisation of the spatial variability (i.e. representing the ‘actual’ in-situ variability). The interval 248
distance between the CPTs in the horizontal direction is Δy = 5 m, and the first CPT is located at y = 249
2.5 m. 250
Figure 2 shows an example realisation, to illustrate the stages involved in constructing the conditional 251
random field. It shows (a) the unconditional field generated using LAS, (b) the Kriged field based on 252
the unconditionally simulated cell values at the measurement locations, (c) the Kriged field based on 253
the measured data (taken from the reference field), and (d) the conditional random field. It can be seen 254
that the conditional field eliminates unrealistic values from the unconditional simulation by honouring 255
the measurement data at the measurement locations (e.g. corresponding to the centre of the dashed 256
circle in the case of the first CPT). The cross-section from which the CPTs were taken is also shown in 257
Figures 2(e) and 2(f), together with the known CPT profiles. It is seen that the known CPT profiles are 258
honoured in the conditional random field. Note that, in order to better visualise the fields, a local 259
colour scale is used for all sub-figures in Figure 2. 260
In order to validate the consistency of the conditioning, the following estimator of the correlation 261
structure along the vertical or horizontal directions of the random field is used to back-figure the 262 covariance structure: 263
(
)
(
)
11
ˆ
(
)
n jˆ
ˆ
j i Z i j Z iC
j
Z
Z
n
j
τ
τ
−µ
+µ
== ∆ =
−
−
−
∑
(11) 264where j = 0, 1, …, n-1, n is the number of data points in the vertical or horizontal direction,
τ
jis the 265lag distance between xi and
x
i j+ , ∆ is the distance between two adjacent cells vertically or τ 266horizontally,
µ
ˆZis the estimated mean, Z is the random soil property and Zi is the sample of Z. The 267correlation function is then ρ τˆ ( )j =Cˆ( )τj Cˆ(0), where Cˆ(0)=σˆ2Zand σˆ2Zis the estimated variance 268
[3]. 269
Figure 3 shows the back-figured (a) vertical and (b) horizontal covariances for the unconditional and 270
conditional random fields averaged over 200 realisations, as well as the sample (i.e. CPT) covariances 271
and exact covariances (i.e. equation (2) with only those terms that are associated with the vertical or 272
horizontal direction in 1D). It can be seen that the conditional random field preserves the covariance 273
structure reasonably well in both the horizontal and vertical directions, and that the correlation 274
function fits well the sample correlation for the first quarter of the data points (i.e. n/4) [5, 9]. It is also 275
seen that the covariance for the conditional field lies in between those for the unconditional field and 276
the sampling points. 277
4. Applications
278Two simple examples concerning slope stability are presented in this section, to illustrate how the 279
technique presented in this paper may be used as an aid to geotechnical design. The first involves 280
finding the optimum locations for CPT profiles, in order to minimise the uncertainty in assessing the 281
reliability of a slope. The second involves a cost-effective design with regard to the slope angle when 282
field measurements have already been made (i.e. the positions where the CPT data were taken are 283
already known). 284
Both examples are presented in terms of the uncertainties in the slope response (with respect to factor 285
of safety). The factors of safety are calculated by 3D finite elements using the strength reduction 286
method [47], with the analyses being undertaken within a probabilistic (RFEM) framework; a 287
flowchart for carrying out such a simulation is shown in Figure 4. The undrained clay behaviour has 288
been modelled using a linear elastic, perfectly plastic Tresca soil model. The clay has a unit weight of 289
20 kN/m3, a Young’s modulus of 100 MPa and a Poisson’s ratio of 0.3. With reference to Figures 5 290
and 13, the finite element boundary conditions are: a fixed base, rollers on the back of the domain 291
preventing displacements perpendicular to the back face, and rollers on the two ends of the domain, 292
allowing only settlements and preventing movements in the other two directions (i.e. the out-of-slope-293
face and longitudinal directions). A full explanation of these boundary conditions is given in Spencer 294
[25] and Hicks and Spencer [19]. 295
The random field cell values are mapped onto the 2×2×2 Gauss points in each 20-node finite element, 296
in order to simulate the spatial variability more accurately [19, 48]. Note that the random fields (both 297
conditional and unconditional) have been mapped onto a finite element mesh with an element aspect 298
ratio equal to 2.0 (see Figure 5) to save time for the finite element analyses [25]. A detailed description 299
of how the random field cell values, in this case based on a cell size of 0.25×0.25×0.25 m, are mapped 300
onto the larger non-cubic finite elements is given in Hicks and Spencer [19]. 301
Note that field test (e.g. CPT) data are not directly used in the following examples. That is, the direct 302
measurements from geotechnical tests are typically not directly applicable in a design. Instead, a 303
transformation model is needed to relate the test measurement (e.g. tip resistance from a CPT test) to 304
an appropriate design property (e.g. the undrained shear strength) [49]. The uncertainty involved in the 305
transformation model is not considered in this paper. 306
4.1 Example 1
307
The first example considers a proposed 45˚, 5 m high, 50 m long slope, that is to be cut from a 308
heterogeneous clay deposit characterised by an undrained shear strength with the following statistics: 309
mean, µ = 20 kPa; standard deviation, σ = 4 kPa; vertical scale of fluctuation, θv = 1.0 m; and 310
horizontal scale of fluctuation, θh = 6.0 m. A question arises as to how to design the sampling strategy 311
for the soil deposit. For example, if 5 CPTs are to be conducted in a straight line along the axis of the 312
proposed slope, where is the best location to site the CPTs such that the designed slope will have the 313
smallest uncertainty in the realised factor of safety F? Hence, this example first investigates the 314
influence of the CPT locations on the standard deviation of the realised factor of safety, followed by 315
the influence of CPT intensity. 316
Figure 5 shows a cross-section through the slope, and 10 possible positions to locate the CPTs (i = 0, 317
1, …, 9). Note that the CPTs are taken to be equally spaced (i.e. at 10 m centres) in the third 318
dimension, and that the first and fifth CPTs are located at 5 m and 45 m along the slope axis (see 319
Figure 9(a)). Furthermore, the CPTs are carried out before the slope is excavated, in a block of soil of 320
dimensions 10×50×5 m as indicated in the figure. 321
Both conditional and unconditional RFEM simulations were carried out, using 500 realisations per 322
simulation, to investigate how the structure response (in this case, the realised factor of safety) 323
changes as the conditioning location changes. Figure 6(a) shows that the uncertainty in the realised 324
factor of safety reduces after conditioning, i.e. after making use of the available CPT information 325
about the soil variability, as indicated by the narrower distribution of realised factor of safety for the 326
conditional simulation. In this figure, the reduction in uncertainty is due to CPT data being taken from 327
location i = 5. 328
Figure 6(b) shows the sampling efficiency indices with respect to the different CPT locations, in which 329
the sampling efficiency index is defined as 330 u se i
I
σ
σ
=
(12) 331 14where
σ
uis the standard deviation of the realised factor of safety for the unconditional simulation, and 332i
σ
is the standard deviation of the realised factor of safety for the conditional simulation based on 333column position i. Hence Ise=1 if the simulation is not conditioned. Clearly, there exists an optimum 334
position (in this case, i = 5) to locate the CPTs; i.e. the uncertainty is a minimum if the CPTs are 335
located along the crest of the proposed slope. In contrast, when i = 0 and i = 1, there is little 336
improvement, because the potential failure planes (in the various realisations) generally pass through 337
zones where the shear strength is, at most, only weakly correlated to values at the left-hand boundary 338
(due to θh being only 6 m in this case). It is interesting to note that, although there is not much 339
information included in the slope stability calculation when i = 9, i.e. for the CPTs at the slope toe, the 340
reduction in uncertainty is still noticeable, due to the CPTs being located in the zone where slope 341
failure is likely to initiate. This observation highlights that the location of additional information may 342
matter more than how much additional information there is (e.g. contrast the large difference in the 343
amount of directly utilised data between CPT locations i = 0 and i = 9). 344
However, it should be remembered that Figures 6(a)-6(b) are for the case of ξ = 6 (corresponding to θh 345
= 6 m) and that ξ often takes a larger value in practice. Figures 6(c)-6(f) show that, for ξ = 12 and ξ = 346
24, the reduction in uncertainty relative to the unconditional case is greater. Moreover, improved 347
values of Ise are obtained for CPT locations near the right and left boundaries, due to the higher 348
correlation of soil properties in the horizontal direction. 349
Figure 7 summarises the results as a function of the degree of anisotropy of the heterogeneity
ξ
. It is 350seen that the best locations for carrying out the 5 CPTs are at i = 5, 6 and 7. As the value of
ξ
351increases, the sampling efficiency indices increase due to the decreasing Kriging variance 2 e
σ
, as 352illustrated in Figure 8 for a y–z slice at i = 5 (i.e. corresponding to where the CPTs are located). It is 353
seen that, for larger values of
ξ
, the Kriging variance between CPTs can drop well below the input 354variance of the shear strength (i.e.
σ
e2≤
16 kPa
2). Moreover, carrying out CPTs at some distance to 355the left or right of the slope crest for higher values of
ξ
can have a similar effect to carrying out CPTs 356near the crest for smaller values of
ξ
. For example, Figure 7 shows that the sampling efficiency index 357for
ξ
=24 at i= is approximately the same as that for 2ξ
=12 at i = 5, 6 and 7. 358Note that the same reference 3D random field is used to represent the ‘real’ field situation in 359
conditioning the random fields in each RFEM analysis. The 3D random fields are conditioned before 360
being mapped onto the finite element mesh, so that they are consistent with sampling the ground 361
before the slope is cut. Hence, for i = 6, 7, 8 and 9, although the CPT measurements are directly used 362
for fewer cells in the FE mesh, they nevertheless have an impact on all cell values via the lateral 363
spatial correlation of soil properties in the original ground profile. 364
If a second row of CPT tests (at position j) is to be performed in a second phase of the site 365
investigation (e.g. as illustrated in Figure 9(b)), the above procedure can be repeated by changing j in 366
the range 0–9 to locate the best positions for the new CPTs, assuming that the position of the first set 367
of CPT profiles has been set to i = 5. This is shown in Figure 10 for the case of
ξ
=6. Figure 10(a) 368shows the probability distributions of the realised factor of safety for the unconditional simulation, the 369
conditional simulation for one row of CPTs at i = 5 and the conditional simulation for an additional 370
row of CPTs at position j = 0. It is seen that the confidence level in the project has been further 371
increased by the second phase of site investigation. Figure 10(b) shows the sampling efficiency indices 372
for various locations j of the second row of CPTs. It suggests that the best location for carrying out the 373
second phase of site investigation can be at either side of the slope crest (at a distance of 374
approximately 3 m (i.e. θh/2) from the crest). 375
To further investigate the influence of CPT intensity on the uncertainty in the realised factor of safety, 376
conditional simulations involving different numbers of CPTs (and thereby different distances (Δ) 377
between adjacent CPTs) have been carried out for the case of ξ = 6, 12 and 24. Figure 11 shows the 378
plan views of CPT layouts for ncpt = 3, 5, 9, 17 and 25 (corresponding to CPT spacings of Δ = 20, 10, 5, 379
3 and 2 m, respectively), with the locations of the CPTs in the x-direction being fixed at i = 5. Figure 380
12 shows the influence of CPT intensity on the sampling efficiency indices for the three values of
ξ
. 381It is seen that there is only a marginal benefit in increasing the scope of the investigation by having 382
CPT spacings less than Δ ≈ θh/2, especially for the
ξ
=6 andξ
=12cases. Forξ
=24, the sampling 383efficiency index is as high as 4 when Δ ≈ θh/2, although more CPTs (i.e. Δ ≈ θh/4, ncpt = 9) may 384
improve the sampling efficiency to a value of 4.5. However, the general finding from Figures 10(b) 385
and 12 is that the optimal sampling distance is around θh/2 for the problem investigated, based on the 386
assumed correlation function. 387
4.2 Example 2
388
In the second example, a soil deposit characterised by spatially varying undrained shear strength is to 389
be excavated to form a slope of a certain angle. Site investigations have been conducted based on CPT 390
tests. The question is: In order to satisfy a target reliability level of, for example, 95%, as suggested in 391
Eurocode [50] and discussed in Hicks and Nuttall [51], how steep should the slope be designed? 392
Figure 13 shows three possible slope angles, with the corresponding finite element mesh 393
discretisations. The slope is 5 m high and 50 m long in the third dimension, and the left-hand boundary 394
is taken to be 15 m from the slope toe. Five CPTs were taken along the length of the slope at 10 m 395
centres, at the location of the column of Gauss points nearest the slope crest for the 1:1 slope, as seen 396
in the figure. The clay soil has a mean undrained shear strength of 21 kPa, a coefficient of variation of 397
0.2, a vertical scale of fluctuation of 1 m and a horizontal scale of fluctuation of 12 m. 398
The three candidate slopes are (vertical:horizontal) 1:2, 1:1 and 2:1. Based on only the mean undrained 399
shear strength, these three slopes have deterministic factors of safety Fd of 1.73, 1.29 and 1.07. Both 400
conditional and unconditional simulations were carried out to investigate the reliability of each slope, 401
and, for each simulation, 500 realisations were analysed. Note that, as in the previous example, one 402
reference random field is generated first and assumed to represent the real field situation. The 403
conditional random fields used in the RFEM analyses are therefore based on CPT measurements taken 404
from this ‘real’ field. 405
The stability of the slopes was calculated by the strength reduction method by applying gravitational 406
loading. The probability density functions of the realised factor of safety are shown in Figure 14 for 407
the three slopes, for both conditional and unconditional simulations. The deterministic factors of safety 408
Fd, i.e. the factors of safety based on the mean property values, are also shown. It is seen that, if 409
unconditional simulation is used, there is a significant chance that the 2:1 slope will fail (the 410
probability of failure is the area under the pdf for the realised factor of safety smaller than 1.0). 411
Unsurprisingly, the gentlest (i.e. 1:2) slope has the lowest probability of failure. However, once again, 412
conditional simulations significantly reduce the uncertainty in the structural response, as clearly 413
demonstrated by the narrower probability distributions. In particular, the reliability of the steepest 414
slope increases from 77% to 99% when the CPT measurements are taken into account. 415
The results show that, if unconditional simulations are used, the 1:1 and 1:2 slopes satisfy a target 416
reliability level of 95%, whereas the 2:1 slope does not. However, when the additional information 417
from the CPT profiles is used, all three cases meet the target reliability. This means that the 418
embankment may be designed to a slope angle of 2:1 if the CPT measurements are used in the 419
simulation, which is, if possible, a more logical thing to do. This has implications for the soil volume 420
to be excavated and thereby cost, although the cost can be site and situation dependent (e.g. on 421
whether there are nearby structures). A best design is a design that meets the requirements set by 422
standards, while, at the same time, minimising the cost. In this case, the steepest slope is likely to be 423
the most cost-effective design. 424
5. Conclusions
425An approach for conditioning 3D random fields based on CPT measurements has been implemented 426
and validated, and then applied to two numerical examples to illustrate its potential use for 427
geotechnical site exploration and cost-effective design. It has been shown that conditional simulations 428
based on CPT data are able to increase the confidence in a design’s success or failure. Indeed, the 429
reliability from a conditional simulation can be thought of as a conditional reliability (or conditional 430
probability of failure not occurring), i.e. based on a ‘posterior’ distribution of the structure 431
performance after taking account of the spatial distribution of all the measured CPT data points. In 432
contrast, the unconditional simulation based on random field theory only results in a ‘prior’ 433
distribution of the structure response. This was clearly demonstrated by the updating of the probability 434
density distributions in the two numerical examples. Although Bayesian updating is not used in this 435
paper, the effect is similar. 436
If further CPT measurements are required, the approach can be repeated for updating the response 437
probability density function. In this way, the confidence in the probability of failure or survival will be 438
further increased. In fact, in many cases a multi-stage site investigation may be carried out, with the 439
results of the initial analysis guiding further field tests. As demonstrated in the first example, if a 440
second stage of site exploration were to be conducted, it is possible to find out the optimum location 441
for the additional testing. This highlights the method’s potential use in directing site exploration 442
programmes and thereby improving the efficient use of field measurements. For the first example 443
considered in this paper, an optimal sampling distance of half the horizontal scale of fluctuation was 444
identified when an exponential correlation function is used. For the second example, the conditional 445
simulation led to a more cost-effective design. 446
Acknowledgements
447This research was funded by the China Scholarship Council (CSC) and by the Section of Geo-448
Engineering at Delft University of Technology. It was carried out on the Dutch National e-449
infrastructure with the support of the SURF Foundation. Special thanks are given to SURFsara advisor 450
Anatoli Danezi for her kind support in developing a computing strategy. 451
References
452[1] DeGroot DJ, Baecher GB. Estimating autocovariance of in-situ soil properties. ASCE Journal of
453
Geotechnical Engineering 1993; 119 (1): 147–66.
454
[2] Fenton GA. Estimation for stochastic soil models. ASCE J. Geotech. Geoenv. Eng. 1999; 125(6): 470–85.
455
[3] Fenton GA. Random field modeling of CPT data. ASCE J. Geotech. Geoenv. Eng. 1999; 125(6): 486–98.
456
[4] Hicks MA, Onisiphorou C. Stochastic evaluation of static liquefaction in a predominantly dilative sand fill.
457
Géotechnique 2005; 55(2): 123–33.
458
[5] Jaksa MB, Kaggwa WS, Brooker PI. Experimental evaluation of the scale of fluctuation of a stiff clay. In:
459
Proc. 8th Int. Conf. on Applications of Statistics and Probability in Civil Engineering, Sydney; 1999. p. 415–22.
460
[6] Lloret-Cabot M, Fenton GA, Hicks MA. On the estimation of scale of fluctuation in geostatistics. Georisk:
461
Assessment and Management of Risk for Engineered Systems and Geohazards 2014; 8(2): 129–40.
462
[7] Lundberg AB, Li Y. Probabilistic characterization of a soft Scandinavian clay supporting a light quay
463
structure. In: Proc. 5th International Symposium on Geotechnical Safety and Risk, Rotterdam; 2015. p.170–75.
464
[8] Phoon KK, Kulhawy FH. Characterization of geotechnical variability. Canadian Geotechnical Journal 1999;
465
36 (4): 612–24.
466
[9] Uzielli M, Vannucchi G, Phoon KK. Random field characterisation of stress-normalised cone penetration
467
testing parameters. Géotechnique 2005; 55(1): 3–20.
468
[10] Zhao HF, Zhang LM, Xu Y, Chang D.S. Variability of geotechnical properties of a fresh landslide soil
469
deposit. Engineering Geology 2013; 166: 1–10.
470
[11] Fenton GA, Griffiths DV. Three-dimensional probabilistic foundation settlement. ASCE J. Geotech.
471
Geoenv. Eng. 2005; 131(2): 232–39.
472
[12] Jaksa MB, Goldsworthy JS, Fenton GA, Kaggwa WS, Griffiths DV, Kuo YL, Poulos HG. Towards reliable
473
and effective site investigations. Géotechnique 2005; 55(2): 109–21.
474
[13] Kuo YL, Jaksa MB, Kaggwa WS, Fenton GA, Griffiths DV, Goldsworthy JS. Probabilistic analysis of
475
multi-layered soil effects on shallow foundation settlement. In: Proc. 9th Australia New Zealand Conference on
476
Geomechanics, Aukland, New Zealand; 2004. p. 541–47.
477
[14] Griffiths DV, Fenton GA. Observations on two- and three-dimensional seepage through a spatially random
478
soil. In: Proc. 7th Int. Conf. on Applications of Statistics and Probability in Civil Engineering, Paris, France;
479
1995. p. 65–70.
480
[15] Griffiths DV, Fenton GA. Three-dimensional seepage through spatially random soil. ASCE J. Geotech.
481
Geoenv. Eng.1997; 123(2): 153–60.
482
[16] Griffiths DV, Fenton GA. Probabilistic analysis of exit gradients due to steady seepage. ASCE J. Geotech.
483
Geoenv. Eng. 1998; 124(9): 789–97.
484
[17] Popescu R, Prevost JH, Deodatis G. 3D effects in seismic liquefaction of stochastically variable soil
485
deposits. Géotechnique 2005; 55(1): 21–31.
486
[18] Griffiths DV, Huang J, Fenton GA. On the reliability of earth slopes in three dimensions. Proc. R. Soc. A
487
2009; 465(2110): 3145–64.
488
[19] Hicks MA, Spencer WA. Influence of heterogeneity on the reliability and failure of a long 3D slope.
489
Comput Geotech 2010; 37: 948–55.
490
[20] Hicks MA, Chen J, Spencer WA. Influence of spatial variability on 3D slope failures. In: Proc. 6th Int. Conf.
491
Computer Simulation Risk Analysis and Hazard Mitigation, Kefalonia; 2008. p. 335–42.
492
[21] Hicks MA, Nuttall JD, Chen J. Influence of heterogeneity on 3D slope reliability and failure consequence.
493
Computers and Geotechnics 2014; 61: 198–208.
494
[22] Li Y, Hicks MA. Comparative study of embankment reliability in three dimensions, In: Proc. 8th European
495
Conference on Numerical Methods in Geotechnical Engineering (NUMGE), Delft; 2014. p. 467–72.
496
[23] Li Y, Hicks MA, Nuttall JD. Probabilistic analysis of a benchmark problem for slope stability in 3D. In:
497
Proc. 3rd Int. Symp. Computational Geomech, Krakow, Poland; 2013. p. 641–8.
498
[24] Li YJ, Hicks MA, Nuttall, JD. Comparative analyses of slope reliability in 3D.Engineering Geology
499
2015;196: 12–23.
500
[25] Spencer WA. Parallel stochastic and finite element modeling of clay slope stability in 3D. PhD thesis,
501
University of Manchester, UK; 2007.
502
[26] Spencer WA, Hicks MA. 3D stochastic modelling of long soil slopes, In: Proc. 14th Conf. of Assoc. for
503
Computational Mechanics in Engineering, Belfast; 2006. p. 119–22.
504
[27] Spencer WA, Hicks MA. A 3D finite element study of slope reliability. In: Proc. 10th Int. Symp. Num
505
Models in Geomech, Rhodes; 2007. p. 539–43.
506
[28] Fenton GA, Griffiths DV. Risk assessment in geotechnical engineering. New York: John Wiley & Sons;
507
2008.
508
[29] Chiles JP, Delfiner P. Geostatistics: modeling spatial uncertainty. John Wiley & Sons; 2009.
509
[30] Lloret-Cabot M, Hicks MA, van den Eijnden AP. Investigation of the reduction in uncertainty due to soil
510
variability when conditioning a random field using Kriging. Géotechnique letters 2012; 2: 123–7.
511
[31] van den Eijnden AP, Hicks MA. Conditional simulation for characterizing the spatial variability of sand
512
state. In: Proc. 2nd International Symposium on Computational Geomechanics, Dubrovnik, Rhodes, Greece;
513
2011. p. 288–96.
514
[32] Vanmarcke EH, Fenton GA. Conditioned simulation of local fields of earthquake ground motion. Structural
515
Safety 1991; 10(1): 247–64.
516
[33] Journel, AG. Geostatistics for conditional simulation of ore bodies. Economic Geology 1974;69(5): 673–87.
517
[34] Olea, RA. Systematic sampling of spatial functions. Series of Spatial Analysis No. 7, Kansas Geological
518
Survey, Lawrence, Kans; 1984.
519
[35] Frimpong, S, Achireko PK. Conditional LAS stochastic simulation of regionalized variables in random
520
fields. Computational Geosciences 1998; 2: 37–45.
521
[36] Journel AG, Huijbregts CJ. Mining Geostatistics. New York: Academic Press; 1978.
522
[37] Fenton GA. Error evaluation of three random field generators. ASCE J. Eng. Mech. 1994; 120 (12): 2478–
523
97.
524
[38] Ji J, Liao HJ, Low BK. Modeling 2D spatial variation in slope reliability analysis using interpolated
525
autocorrelations. Computers and Geotechnics 2012; 40: 135–46.
526
[39] Phoon KK, Huang SP, Quek ST. Simulation of second-order processes using Karhunen–Loeve expansion.
527
Computers & Structures 2002; 80(12): 1049–60.
528
[40] Fenton GA, Vanmarcke EH. Simulation of random fields via local average subdivision. ASCE Journal of
529
Engineering Mechanics 1990; 116(8): 1733–49.
530
[41] Krige DG. A statistical approach to some basic mine valuation problems on the Witwatersrand. Journal of
531
the Chemical, Metallurgical and Mining Society of South Africa 1951; 52 (6): 119–39.
532
[42] Cressie N. The origins of Kriging. Mathematical Geology 1990; 22 (3): 239–52.
533
[43] Wackernagel H. Multivariate geostatistics: An introduction with applications. Germany: Springer; 2003.
534
[44] Fenton GA. Simulation and analysis of random fields. PhD Thesis, Princeton University; 1990.
535
[45] Vanmarcke EH. Random fields: analysis and synthesis. Cambridge, Massachusetts: The MIT Press; 1983.
536
[46] Fenton GA. Data analysis/geostatistics. In: Probabilistic methods in geotechnical engineering. Griffiths DV,
537
Fenton GA, editors. New York: Springer; 2007. p. 51–73.
538
[47] Smith IM, Griffiths DV, Margetts L. Programming the finite element method, 5th ed. New York: John
539
Wiley & Sons; 2013.
540
[48] Hicks MA, Samy K. Influence of heterogeneity on undrained clay slope stability. Quarterly Journal of
541
Engineering Geology and Hydrogeology 2002; 35: 41–9.
542
[49] Phoon KK, Kulhawy FH. Evaluation of geotechnical property variability. Canadian Geotechnical Journal
543
1999; 36(4): 625–39.
544
[50] European Committee for Standardisation (CEN). Eurocode 7: geotechnical design. Part 1: general rules. EN
545
1997-1, CEN; 2004.
546
[51] Hicks MA, Nuttall JD. Influence of soil heterogeneity on geotechnical performance and uncertainty: a
547
stochastic view on EC7. In: Proc. 10th International Probabilistic Workshop, Stuttgart; 2012. p. 215–27.
548 549
550 551 552 553 554
Appendix
555A.1 Forming the left-hand-side matrix for Kriging
556
Suppose there are k×m CPT locations that follow a rectangular grid at the ground surface. That is, 557
there are k rows in the x direction and, within each row, m CPT profiles in the y direction (Figure 1). 558
Assuming that there are n data points for each CPT profile, the global numbering scheme for all the 559
CPT data points is shown in Figure A.1 for the case of k = 2. 560
Following the basic equation (6), of size N + 1 = k×m×n + 1, the left-hand-side matrix is formulated 561 as 562 1,1 1,2 1,3 1, 1, 1 1, 2 1, 3 1,2 2,1 2,2 2,3 2, 2, 1 2, 2 2, 3 2,2 3,1 3,2 3,3 3, 3, 1 3, 2 3, 3 3,2 ,1 ,2 ,3 , , 1 , 2 , 3 ,2
1
1
1
1
1
1
1
m m m m m m m m m m m m m m m lhs km km km km m km m km m km m km m + + + + + + + + + + + +=
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
γ
v
v
v
v
v
v
v
v
2
2
2
2
2
2
2
2
2
2
2
2
2
1,( 1) 1 1,( 1) 2 1,( 1) 3 1, 2,( 1) 1 2,( 1) 2 2,( 1) 3 2, 3,( 1) 1 3,( 1) 2 3,( 1) 3 3, ,( 1) 1 ,( 1) 2 ,( 1) 3 ,1
1
1
1
1
1
1
1
1
0
k m k m k m km k m k m k m km k m k m k m km km k m km k m km k m km km − + − + − + − + − + − + − + − + − + − + − + − +
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
2
2
2
2
2
2
2
(A1) 563in which
v
i j, is a matrix representing the correlation structure between CPTi and CPTj (where each 564CPT has n data points), 565
( 1) 1,( 1) 1 ( 1) 1,( 1) 2 ( 1) 1,( 1) 3 ( 1) 1,( 1) ( 1) 2,( 1) 1 ( 1) 2,( 1) 2 ( 1) 2,( 1) 3 ( 1) 2,( 1) ( 1) 3,( 1) 1 ( 1) 3,( 1) 2 ( 1) 3,( , i n j n i n j n i n j n i n j n n i n j n i n j n i n j n i n j n n i n j n i n j n i n j i j
d
d
d
d
d
d
d
d
d
d
d
− + − + − + − + − + − + − + − + − + − + − + − + − + − + − + − + − + − + − + − + − + −=
v
2
2
1) 3 ( 1) 3,( 1) ( 1) ,( 1) 1 ( 1) ,( 1) 2 ( 1) ,( 1) 3 ( 1) ,( 1) n i n j n n i n n j n i n n j n i n n j n i n n j n nd
d
d
d
d
+ − + − + − + − + − + − + − + − + − + − +
2
2
(A2) 566where (i, j) = 1, 2, 3, …, m, m+1, m+2, m+3, …, 2m, ……, (k-1)m+1, (k-1)m+2, (k-1)m+3, …, km and 567
dr,s (r = (i-1)n+1, …, (i-1)n+n) (s = (j-1)n+1, …, (j-1)n+n) are the components of the submatrix
v
i j, , 568which can be expressed in the form of a covariance function between data points r and s (equation (2)). 569
A.2 Forming the right-hand-side vector for Kriging
570
The right-hand-side vector is formulated as 571 1 2 3
1
rhs km
=
v
v
v
γ
v
(A3) 572in which
v
pis a vector representing the correlation structure between the estimation point and CPTp, 573 ( 1) 1 ( 1) 2 ( 1) 3 ( 1) p n p n p n p p n nd
d
d
d
− + − + − + − +
=
v
(A4) 574 where p = 1, 2, 3, …, m, m+1, m+2, m+3, …, 2m, ……, (k-1)m+1, (k-1)m+2, (k-1)m+3, …, km and dt 575(t = (p-1)n+1, …, (p-1)n+n) are the components of the subvector
v
p, which can be expressed in the 576form of a covariance function (equation (2)) between data points t and the point at which the value is 577
to be estimated (Figure A.1). 578
The unknown weight vector is 579
1 2 3 x km
µ
=
λ
λ
λ
λ
λ
(A5) 580in which
λ
qis the weight subvector for CPTq, 581 ( 1) 1 ( 1) 2 ( 1) 3 ( 1) q n q n q n q q n nλ
λ
λ
λ
− + − + − + − +
=
λ
(A6) 582 where q = 1, 2, 3, …, m, m+1, m+2, m+3, …, 2m, ……, (k-1)m+1, (k-1)m+2, (k-1)m+3, …, km. 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 25List of Figures
599600
Figure 1. Example CPT sampling strategy (k = 2, m = 5) 601
Figure 2. Example illustrations of the unconditional random field (a), the Kriged field based on the 602
randomly simulated data (b), the Kriged field based on the CPT data (c), the conditional random field 603
(d), cross-sections (e and f) in the longitudinal direction taken from the Kriged field (c) and from the 604
conditional random field (d), respectively. Dashed circle indicates the position of the first CPT in 605
subfigures (a) and (c-d) 606
Figure 3. Vertical and horizontal covariance functions averaged over 200 realisations (θv = 1.0 m, θh = 607
3.0 m) 608
Figure 4. Flowchart for conditional RFEM simulation 609
Figure 5. Finite element mesh and possible numbered CPT locations at a cross-section through the 610
proposed 50 m long slope (dashed lines indicate the excavated soil mass and numbers correspond to 611
Gauss point locations within the finite elements) 612
Figure 6. Simulation results for Example 1 (based on θv = 1.0 m and 500 realisations per simulation) 613
Figure 7. Sampling efficiency indices for various values of
ξ
614Figure 8. Kriging variance for various values of
ξ
(y–z slice at i = 5) 615Figure 9. CPT layout illustration (plan view) for a single row (a) and two rows (b) 616
Figure 10. Influence of CPT location j during the second phase of site investigation (based on θv = 1.0 617
m and 500 realisations per simulation) 618
Figure 11. CPT layouts (plan views) for various numbers of boreholes (ncpt = 3, 5, 9, 17, 25 and Δ 619
denotes the distance between CPTs) 620
