# CPT-Based Probabilistic Characterization of Undrained Shear Strength of Clay

Authors: Mi Tian; Xiaotao Sheng
Journal: Advances in Civil Engineering
 (2020)

---

## Abstract

Applying random field theory involves two important issues: the statistical homogeneity (or stationarity) and determination of random field parameters and correlation function. However, the profiles of soil properties are typically assumed to be statistically homogeneous or stationary without rigorous statistical verification. It is also a challenging task to simultaneously determine random field parameters and the correlation function due to a limited amount of direct test data and various uncertainties (e.g., transformation uncertainties) arising during site investigation. This paper presents Bayesian approaches for probabilistic characterization of undrained shear strength using cone penetration test (CPT) data and prior information. Homogeneous soil units are first identified using CPT data and subsequently assessed for weak stationarity by the modified Bartlett test to reject the null hypothesis of stationarity. Then, Bayesian approaches are developed to determine the random field parameters and simultaneously select the most probable correlation function among a pool of candidate correlation functions within the identified statistically homogeneous layers. The proposed approaches are illustrated using CPT data at a clay site in Shanghai, China. It is shown that Bayesian approaches provide a rational tool for proper determination of random field model for probabilistic characterization of undrained shear strength with consideration of transformation uncertainty.


---

## Body

## 1. Introduction

Soil materials are nature materials, and their properties inherently vary from one location to another due to their natural geologic process, which is known as “inherent spatial variability (ISV).” ISV of design soil properties (e.g., undrained shear strength,Su) can be probabilistically characterized using random field theory in probability-based geotechnical analyses and designs. In the context of random field theory, a design soil property within a statistically homogeneous layer is described by a series of random variables. These random variables are with the same mean, μ, and standard deviation, σ; and the autocorrelation among these random variables depends on scale of fluctuation, λ, and the correlation function (such as those shown in Figure 1) (e.g., [2–6]). For a statistically homogeneous soil layer, the random field model of design soil properties can be uniquely defined by a given set of random field parameters (i.e., μ, σ, and λ) and correlation function. Note that applying random field model to describe ISV of design soil properties involves two important issues: the statistical homogeneity (or stationarity) and determination of random field parameters and correlation function.Figure 1
Four commonly used correlation functions in geotechnical engineering (after [1]).Statistical homogeneity (or stationarity) means that the mean and standard deviation of design soil properties are assumed to be spatially constant, and the correlation function only depends on their separation distance between two different locations (e.g., [2, 7, 8]). Statistical homogeneity is an important prerequisite for statistical analyses (e.g., determining random field parameters) on a set of observation data (e.g., [2, 7, 9, 10]). If the stationarity of observation data is not properly verified, the random field model may be biased estimate. However, the profiles of soil properties are typically assumed to be statistically homogeneous or stationary without rigorous statistical verification in the previous studies (e.g., [11, 12]).For observation data sets that satisfy the statistically homogeneous or stationary conditions, the inherent soil variability can be evaluated using random field model. However, it is a challenging task to identify random field parameters and select the proper correlation function of design soil properties at a given site, because only a small portion of soil materials are examined during site investigation and the amount of observation data is usually too limited for meaningful estimates of random field parameters. Furthermore, design soil properties (e.g., undrained shear strength,Su) used in geotechnical designs are often estimated indirectly through transformation models (e.g., the empirical regression between the cone penetration test (CPT) data and undrained shear strength, Su). The transformation model is not a perfect relationship but is associated with uncertainties/dispersion about a mean trend, namely, “transformation uncertainty” (e.g., [13]), which shall be rationally considered to describe ISV of design soil properties. This can be addressed as an inverse analysis problem within a Bayesian framework, which can provide a rational vehicle to determine the random field parameters and correlation function by combining the project-specific data with prior information (e.g., [14]). Bayesian approaches have been increasingly used for probabilistic characterization of soil properties, such as probabilistic characterization of sand effective friction angle (e.g., [11, 12, 15]) and normalized cone tip resistance (e.g., [16]).In this study, Bayesian approaches are adapted to characterize the ISV of undrained shear strength of clay using CPT data. The homogeneous soil units are first identified based on the site-specific CPT data and subsequently assessed for stationarity using the modified Bartlett test. Then, Bayesian approaches are developed to determine the random field parameters (i.e.,μ, σ, and λ) and simultaneously select the most probable correlation function among the candidate correlation functions (as shown in Figure 1) of undrained shear strength within a statistically homogeneous layer. Prior knowledge and project-specific CPT data are incorporated into the proposed Bayesian approaches to describe the ISV of undrained shear strength. Finally, the Bayesian approaches are illustrated using a set of real-life CPT data.


## 2. Identification of Statistical Homogeneity or Stationarity

### 2.1. Identification of Homogeneous Soil Units

An important prerequisite for the determination of random field model is the condition of observation data, if necessary, to achieve statistical homogeneity (or stationarity). The mean and standard deviation tend to be specifically related to the soil type, especially for the mean. The correlation function of soil properties also depends on the soil type (e.g., [17]). Investigating uniform soil units in composition or behavior is of great significance for the assumption of statistical homogeneity, which should be satisfied in the subsequent statistical analyses (e.g., [1, 7, 10, 17–21]. Homogeneous soil units can be examined using CPT data, which provides direct information regarding the soil behavior type by measuring both the cone tip resistance qc and sleeve friction fs. Based on the measured qc and fs, the soil behavior type can be identified using the soil behavior classification index, Ic (e.g., [19, 20, 22]. Consider, for example, the soil behavior classification index shown in Table 1. The soil behavior type is divided into nine different types corresponding to the ranges of Ic. Based on which zones the measured Ic is located in, the corresponding soil type is determined accordingly. Ic is calculated by the normalized cone resistance, Qt, and the normalized friction ratio, FR:(1)Ic=3.47−logQt2+logFR+1.2220.5.Table 1
Soil behavior type indexIc (after [20]).
Soil behavior zone
I
c range
Soil behavior type
1
Nonappl.
Sensitive, fine grained
2
I
c > 3.60
Organic soils—peats
3
2.95 <Ic < 3.60
Clay, silty clay
4
2.60 <Ic < 2.95
Silt mixtures—clayey silt to silty clay
5
2.05 <Ic < 2.60
Sand mixtures—silty sand to sandy silt
6
1.31 <Ic < 2.05
Clean sand, silty sand
7
I
c < 1.31
Gravelly sand to sand
8
Nonappl.
Very stiff sand to clayey sand
9
Nonappl.
Very stiff, fine grainedThe normalized cone resistanceQt is proposed by Robertson [23] and Robertson [24] and it is written as(2)Qt=qt−σv0σv0′,where qt=qc+1−au is the corrected cone tip resistance, in which u is pore pressure measured between the cone tip and the friction sleeve, and a is net area ratio; σv0 is the total vertical stress; σv0′ is the vertical effective stress. The normalized friction ratio FR is given by (e.g., [23, 24])(3)FR=100fsqt−σv0.Homogeneous soil units are identified from the CPT data by delineating the data into sections where coefficient of variation (COV) ofIc is less than 0.1 (e.g., [19, 25]) and Ic does not exhibit gradual changes with depth. These conditions are imposed to ensure that there are sufficient data points for stationarity assessment and the sampling size is large enough for meaningful statistical analyses.


### 2.2. Assessment of Stationarity

The homogeneous soil units are further screened to reject the null hypothesis of stationarity. Due to the limitation of sampling size in geotechnical practice, stationarity generally can only be assessed in a weak sense. Numerous classical statistical tests have been applied to assess weak stationarity, such as statistical runs test, Kendall’s tau test, and Spearman’s rank coefficient, to correlated data (e.g., [20, 26]), which may result in biased assessments. The modified Bartlett statistic (MBS) procedure proposed by Phoon et al. [1] provides a more rational basis for rejecting the null hypothesis of weak stationarity in the correlated case. In the MBS procedure, a continuous Bartlett statistic profile, Bstat, is first generated by moving a sampling window over the soil profile. The peak value of the Bartlett statistic profile, Bmax, is compared with the critical value, Bcrit, for rejection or acceptance of the weak stationarity hypothesis. Phoon et al. [1] has proposed the rejection criterion corresponding to the commonly used correlation functions (shown in Table 2) (e.g., [1]).Table 2
Critical modified Bartlett test statistic at 5% level of significance for various common autocorrelation models (I2 = 1) (after [1]).
Correlation function
Rejection criteria,Bcrit
SECF
0.23
υ
+
0.71
ln
I
1
+
0.91
υ
+
0.23
BNCF
0.30
υ
+
0.29
ln
I
1
+
1.151
υ
−
0.52
SMCF
0.42
υ
−
0.07
ln
I
1
+
2.041
υ
−
3.32
SQECF
0.73
υ
−
0.98
ln
I
1
+
2.351
υ
−
2.45
Note:I1 = Nn/υ is the normalized sampling length, where Nn is the total number of points in entire soil record, υ = λ/∆ is the number of points in one scale of fluctuation, and ∆ is the sampling interval; I2=w/υ is the normalized segment length.Bartlett’s test is used to test the equality of multiple sample variances for independent data sets. For the case of two sample variances, the sampling window is divided into two equal segmentss12 and s22. The Bartlett test statistic is written as(4)Bstat=2.30259w−1A2logs2−logs12+logs22,in which w is the number of data points used to evaluate s12 or s22; s2=s12+s22/2 is the total variance; A=1+1/2w−2 is a constant. The Bartlett statistic profile can be generated by moving a sampling window over the soil profiles. Then, the peak value Bmax is compared with the critical value, Bcrit, to assess the stationarity. The null hypothesis of stationarity is rejected at the 5% level of significance if Bmax > Bcrit.


## 2.1. Identification of Homogeneous Soil Units

An important prerequisite for the determination of random field model is the condition of observation data, if necessary, to achieve statistical homogeneity (or stationarity). The mean and standard deviation tend to be specifically related to the soil type, especially for the mean. The correlation function of soil properties also depends on the soil type (e.g., [17]). Investigating uniform soil units in composition or behavior is of great significance for the assumption of statistical homogeneity, which should be satisfied in the subsequent statistical analyses (e.g., [1, 7, 10, 17–21]. Homogeneous soil units can be examined using CPT data, which provides direct information regarding the soil behavior type by measuring both the cone tip resistance qc and sleeve friction fs. Based on the measured qc and fs, the soil behavior type can be identified using the soil behavior classification index, Ic (e.g., [19, 20, 22]. Consider, for example, the soil behavior classification index shown in Table 1. The soil behavior type is divided into nine different types corresponding to the ranges of Ic. Based on which zones the measured Ic is located in, the corresponding soil type is determined accordingly. Ic is calculated by the normalized cone resistance, Qt, and the normalized friction ratio, FR:(1)Ic=3.47−logQt2+logFR+1.2220.5.Table 1
Soil behavior type indexIc (after [20]).
Soil behavior zone
I
c range
Soil behavior type
1
Nonappl.
Sensitive, fine grained
2
I
c > 3.60
Organic soils—peats
3
2.95 <Ic < 3.60
Clay, silty clay
4
2.60 <Ic < 2.95
Silt mixtures—clayey silt to silty clay
5
2.05 <Ic < 2.60
Sand mixtures—silty sand to sandy silt
6
1.31 <Ic < 2.05
Clean sand, silty sand
7
I
c < 1.31
Gravelly sand to sand
8
Nonappl.
Very stiff sand to clayey sand
9
Nonappl.
Very stiff, fine grainedThe normalized cone resistanceQt is proposed by Robertson [23] and Robertson [24] and it is written as(2)Qt=qt−σv0σv0′,where qt=qc+1−au is the corrected cone tip resistance, in which u is pore pressure measured between the cone tip and the friction sleeve, and a is net area ratio; σv0 is the total vertical stress; σv0′ is the vertical effective stress. The normalized friction ratio FR is given by (e.g., [23, 24])(3)FR=100fsqt−σv0.Homogeneous soil units are identified from the CPT data by delineating the data into sections where coefficient of variation (COV) ofIc is less than 0.1 (e.g., [19, 25]) and Ic does not exhibit gradual changes with depth. These conditions are imposed to ensure that there are sufficient data points for stationarity assessment and the sampling size is large enough for meaningful statistical analyses.


## 2.2. Assessment of Stationarity

The homogeneous soil units are further screened to reject the null hypothesis of stationarity. Due to the limitation of sampling size in geotechnical practice, stationarity generally can only be assessed in a weak sense. Numerous classical statistical tests have been applied to assess weak stationarity, such as statistical runs test, Kendall’s tau test, and Spearman’s rank coefficient, to correlated data (e.g., [20, 26]), which may result in biased assessments. The modified Bartlett statistic (MBS) procedure proposed by Phoon et al. [1] provides a more rational basis for rejecting the null hypothesis of weak stationarity in the correlated case. In the MBS procedure, a continuous Bartlett statistic profile, Bstat, is first generated by moving a sampling window over the soil profile. The peak value of the Bartlett statistic profile, Bmax, is compared with the critical value, Bcrit, for rejection or acceptance of the weak stationarity hypothesis. Phoon et al. [1] has proposed the rejection criterion corresponding to the commonly used correlation functions (shown in Table 2) (e.g., [1]).Table 2
Critical modified Bartlett test statistic at 5% level of significance for various common autocorrelation models (I2 = 1) (after [1]).
Correlation function
Rejection criteria,Bcrit
SECF
0.23
υ
+
0.71
ln
I
1
+
0.91
υ
+
0.23
BNCF
0.30
υ
+
0.29
ln
I
1
+
1.151
υ
−
0.52
SMCF
0.42
υ
−
0.07
ln
I
1
+
2.041
υ
−
3.32
SQECF
0.73
υ
−
0.98
ln
I
1
+
2.351
υ
−
2.45
Note:I1 = Nn/υ is the normalized sampling length, where Nn is the total number of points in entire soil record, υ = λ/∆ is the number of points in one scale of fluctuation, and ∆ is the sampling interval; I2=w/υ is the normalized segment length.Bartlett’s test is used to test the equality of multiple sample variances for independent data sets. For the case of two sample variances, the sampling window is divided into two equal segmentss12 and s22. The Bartlett test statistic is written as(4)Bstat=2.30259w−1A2logs2−logs12+logs22,in which w is the number of data points used to evaluate s12 or s22; s2=s12+s22/2 is the total variance; A=1+1/2w−2 is a constant. The Bartlett statistic profile can be generated by moving a sampling window over the soil profiles. Then, the peak value Bmax is compared with the critical value, Bcrit, to assess the stationarity. The null hypothesis of stationarity is rejected at the 5% level of significance if Bmax > Bcrit.


## 3. Random Field Model of Undrained Shear Strength

Random field model (e.g., [3]) is applied in this paper to describe the ISV of undrained shear strength, Su, within a single statistically homogeneous layer. Insights from soil mechanics show that Su itself is not a fundamental soil property but is related to its corresponding vertical effective stress, σv0′ (e.g., [27–30]). Then, the undrained shear strength ratio r=Su/σv0′, where Su and σv0′ represent the undrained shear strength and vertical effective stress at a given depth, respectively. Normalization by the vertical effective stress is used to remove the trend in the undrained shear strength profile, which generally increases with depth [31]. Other detrending methods can also be applied, such as trend line regression analysis using either linear or polynomial curve fittings [32].Consider, for example,r at different depths are modeled by a series of spatially correlated lognormal variables with a mean μ and standard deviation σ (i.e., one-dimensional stationary lognormal random field) without explicit consideration of the variation trend at different depths. The random field model adopted in this study is generally reasonable for soft clays (e.g., normal consolidated or slightly overconsolidated clays). The logarithm of r (i.e., γ=lnr) is a normal random variable with a mean μN=lnμ−σN2/2 and standard deviation σN=ln1+σ/μ2. The spatial correlation between variations of γ at different depths is specified by the scale of fluctuation λ and a correlation function M. Examples of correlation functions include the single exponential correlation function (SECF), binary noise correlation function (BNCF), second-order Markov correlation function (SMCF), and squared exponential correlation function (SQECF), as shown in Figure 1.Letγ¯=γD1,γD2,…,γDn be a vector of γD at different depths D1, D2, …, Dn. Because γ is considered to be normally distributed with a mean μN and standard deviation σN, γ¯ is a Gaussian vector with a mean vector μNl¯ and covariance matrix σN2R¯, in which l¯ is a vector with n components that are all equal to one and R¯ is the correlation matrix of γD. Then, γ¯ can be written as(5)γ¯=μNl¯+σNL¯TZ¯,in which Z¯ is a standard normal random variable vector; L¯ is an n by n upper-triangular matrix obtained by Cholesky decomposition of R¯. The (i, j)th entry of R¯ represents the correlation coefficient between γDi and γDj at respective depths of Di and Dj, and it is given by the scale of fluctuation λ and a correlation function M as shown in Figure 1.The undrained shear strengthSu of clay can be estimated from in situ test, such as CPT through empirical regression functions. Consider, for example, the following empirical regression function (e.g., [33]):(6)Su=0.0789qt−σv0′.Equation (6) can be rewritten as r=Su/σv0′=0.0789qt−σv0′/σv0′, where the effects of vertical effective stress, σv0′, are taken into account. In a log scale, the logarithm of r is rewritten as(7)lnSuσv0′=lnqt−σv0′σv0′−2.540=lnQt−2.540.A simple Bayesian argument can transform the above equation into(8)lnQt=alnr+b+ε,where a = 1; b = 2.540; and ε is a Gaussian random variable with a mean of zero and a standard deviation σε of 0.34. ε represents the transformation uncertainty associated with the empirical regression function.Letξ¯=ξD1,ξD2,…,ξDnT be a vector of lnQt data at different depths D1, D2, …, Dn. Using equations (5) and (8) gives(9)ξ¯=aμN+bl¯+aσNL¯TZ¯+σεε¯.When the spatial variability is assumed to be independent of transformation uncertainty, that is,Z¯ and ε¯ are independent, ξ¯ is a Gaussian random vector with a mean vector aμN+bl¯ and covariance matrix C¯=a2σN2R¯+σε2I¯, where I¯ is an n × n identity matrix.Knowledge onμ, σ, and λ and correlation function M are required to define the probability distribution of ξ¯. The next section presents Bayesian approaches for identifying the random field parameters (i.e., μ, σ, and λ) and select the most probable correlation function among a pool of candidates (e.g., those shown in Figure 1) based on both project-specific data (i.e., CPT data ξ¯) and prior knowledge.


## 4. Bayesian Identification of Random Field Model

### 4.1. Bayesian Identification of Random Field Parameters

For a given correlation functionM, the stationary normal random field of r within a single statistically homogeneous layer is uniquely represented by the random field parameters X=μ,σ,λ. Note that the correlation function M is assumed to take a specific form (e.g., one of the correlation functions shown in Figure 1) in this section, but its corresponding λ value is unknown herein. A proper form of the correlation function will be determined among a pool of the candidates (e.g., those shown in Figure 1) by a Bayesian model class selection approach in Section 4.2 entitled “Bayesian Selection of Correlation Function.”For a given set of prior knowledge and site observation dataξ¯, there are various possible values of random field parameters, and their respective plausibility can be quantified by the posterior distribution PXξ¯,M under a Bayesian framework, where the condition on M indicates that the correlation function M is assumed to take a specific form. Using Bayes’ theorem, the posterior distribution of random field parameters PXξ¯,M is written as (e.g., [14, 15, 34–36])(10)PXξ¯,M=K−1Pξ¯X,MPXM,in which ξ¯=ξD1,ξD2,…,ξDnT is a set of values of ξD measured at different depths D1, D2, …, Dn in a single clay layer, respectively; K=Pξ¯M is a normalizing constant that is independent of X for a given M, and it is calculated as ∫Pξ¯X,MPXMdX; Pξ¯X,M is the likelihood function that reflects the model fit with observation data ξ¯ for a given M and random field parameters X; and PXM is the prior distribution of random field parameters X reflecting the available prior knowledge on X in the absence of data. When there is no prevailing knowledge on X, a relatively noninformative prior distribution can be used to reflect prior knowledge, for example, a joint uniform distribution that is equal to a constant (e.g., [36]). When local information (e.g., local experience and engineering judgments) is available at a site, the prior knowledge becomes more informative. More sophisticated types of prior distributions (e.g., an arbitrary histogram type of prior distribution, [37]) can be obtained. This may complicate the posterior distribution in equation (10) and lead to difficulties in expressing it analytically and explicitly.As mentioned before,ξ¯ is a Gaussian random vector with a mean vector aμN+bl¯ and covariance matrix C¯. The likelihood function is thus given by(11)Pξ¯X,M=2π−n/2detC¯−1/2exp−12ξ¯−aμN+bl¯TC¯−1ξ¯−aμN+bl¯.The prior distribution and the likelihood function are then substituted into equation (10) to obtain the posterior distribution PXξ¯,M. The posterior distribution involves high dimension integration, which is one key limitation of Bayesian methods. To avoid the computational complexity, Markov Chain Monte Carlo Simulation (MCMCS) is adapted in this study to generate a large number of samples of random field parameters X from the posterior distribution (i.e., equation (10)). Note that MCMCS provides an appropriate tool to obtain the posterior knowledge in Bayesian analysis (e.g., [38]). Metropolis-Hastings (M-H) algorithm is used to generate a large number of X samples to describe the posterior distribution (i.e., equation (10), e.g., [15]). Samples of X are subsequently used to estimate the statistics of random field parameters through conventional statistical analyses.


### 4.2. Bayesian Selection of Correlation Function

The correlation functionM within a single statistically homogeneous clay layer is considered as a prescribed one with unknown λ in the previous section. This section presents a Bayesian model class selection method to determine the most probable correlation function M∗ among a pool of candidate correlation functions for a given set of ξ¯. Consider a number (i.e., NCF) of candidate correlation functions Mk, k = 1, 2, …, NCF within the clay layer, such as SECF, BNCF, SMCF, and SQECF shown in Figure 1. For a given set of ξ¯, the plausibility of the k-th candidate correlation function Mk in the clay layer is quantified by its occurrence probability PMkξ¯. Using Bayes’ theorem, PMkξ¯ is given by (e.g., [15, 36])(12)PMkξ¯=Pξ¯MkPMkPξ¯,k=1,2,…,NCF,where Pξ¯ is a normalizing constant that is independent of Mk and it is calculated as Pξ¯=∑k=1NCFPξ¯MkPMk; PMk is the prior probability of Mk; and Pξ¯Mk is the conditional probability density function of ξ¯ for a given correlation function Mk, and it quantifies the information on Mk provided by CPT data and is referred to as “evidence” for Mk provided by ξ¯.In the case of no prevailing prior knowledge, theNCF candidate correlation functions in the clay layer are considered to have the same prior probability PMk=1/NCF. Since Pξ¯ and PMk in equation (12) are constants, PMkξ¯ is proportional to Pξ¯Mk. The most probable correlation function M∗ is the one that has the maximum value of PMkξ¯, and it also has the maximum value of Pξ¯Mk. Note that Pξ¯Mk is actually the normalizing constant K in equation (10) by setting M as Mk, which involves a multidimensional integration. As discussed in the previous section, there is no need to calculate K (i.e., Pξ¯Mk) when MCMCS is adapted to the posterior distribution in equation (10) for Mk. However, it is needed herein for Bayesian selection of correlation function. MCMCS samples of X generated from equation (10) for Mk and Gaussian copula are used to calculate the evidence Pξ¯Mk (e.g., [15]), maintaining the advantage of no need to directly calculate the multidimensional integration (i.e., K).Using the Theorem of Total Probability, the evidencePξ¯Mk for Mk is expressed as(13)Pξ¯Mk=∫Pξ¯X,MkPXMkdX,where Pξ¯X,Mk is the likelihood function by setting M = Mk; and PXMk is the prior distribution of X. Then, the logarithm of evidence (i.e., lnPξ¯Mk) can be written as (e.g., [36])(14)lnPξ¯Mk=∫lnPξ¯X,MkPXξ¯,MkdX−∫lnPXξ¯,MkPXMkPXξ¯,MkdX.Based on theNp random samples of X estimated from PXξ¯,Mk, lnPξ¯Mk can be estimated as(15)lnPξ¯Mk≈∑i=1NplnPξ¯xi,MkNp−∑i=1NplnPxiξ¯,Mk/PxiMkNp,in which PxiMk is the prior distribution evaluated at i-th sample xi = [μi, σi, λi] of X, which is assumed to be a constant; Pξ¯xi,Mk is the likelihood function calculated using equation (11) for each xi; Pxiξ¯,Mk is posterior distribution evaluated using equation (10) at i-th sample xi.EvaluatingPxiξ¯,Mk needs the posterior distribution PXξ¯,Mk in equation (10), which is unknown and is numerically represented by Np MCMCS samples of X. To evaluate PXξ¯,Mk, copula theory is used to reconstruct PXξ¯,Mk based on Np samples of X in this paper. Details of calculation of the evidence using MCMCS and Gaussian copula can be found in Tian et al. [15].For a givenMk, the MCMCS-based Bayesian approach in the previous subsection is applied to generate random samples of X from the posterior distribution. These samples are subsequently used in equation (15) to evaluate lnPξ¯Mk. The procedure is repetitively performed NCF times for Mk, k = 1, 2, …, NCF to obtain their respective values of lnPξ¯Mk. Then, the most probable correlation function M∗ is determined by comparing the values of lnPξ¯Mk in each candidate correlation function. The one with the maximum value of lnPξ¯Mk is taken as M∗. Its corresponding posterior statistics (e.g., posterior mean values μ′, σ′, and λ′) estimated from MCMCS samples constitute the posterior knowledge on random field parameters. In other words, the most probable correlation function M∗ and posterior estimates (e.g., μ′, σ′, and λ′) of random field parameters are determined simultaneously for probabilistic characterization of the ISV of undrained shear strength based on both prior knowledge and CPT data.


## 4.1. Bayesian Identification of Random Field Parameters

For a given correlation functionM, the stationary normal random field of r within a single statistically homogeneous layer is uniquely represented by the random field parameters X=μ,σ,λ. Note that the correlation function M is assumed to take a specific form (e.g., one of the correlation functions shown in Figure 1) in this section, but its corresponding λ value is unknown herein. A proper form of the correlation function will be determined among a pool of the candidates (e.g., those shown in Figure 1) by a Bayesian model class selection approach in Section 4.2 entitled “Bayesian Selection of Correlation Function.”For a given set of prior knowledge and site observation dataξ¯, there are various possible values of random field parameters, and their respective plausibility can be quantified by the posterior distribution PXξ¯,M under a Bayesian framework, where the condition on M indicates that the correlation function M is assumed to take a specific form. Using Bayes’ theorem, the posterior distribution of random field parameters PXξ¯,M is written as (e.g., [14, 15, 34–36])(10)PXξ¯,M=K−1Pξ¯X,MPXM,in which ξ¯=ξD1,ξD2,…,ξDnT is a set of values of ξD measured at different depths D1, D2, …, Dn in a single clay layer, respectively; K=Pξ¯M is a normalizing constant that is independent of X for a given M, and it is calculated as ∫Pξ¯X,MPXMdX; Pξ¯X,M is the likelihood function that reflects the model fit with observation data ξ¯ for a given M and random field parameters X; and PXM is the prior distribution of random field parameters X reflecting the available prior knowledge on X in the absence of data. When there is no prevailing knowledge on X, a relatively noninformative prior distribution can be used to reflect prior knowledge, for example, a joint uniform distribution that is equal to a constant (e.g., [36]). When local information (e.g., local experience and engineering judgments) is available at a site, the prior knowledge becomes more informative. More sophisticated types of prior distributions (e.g., an arbitrary histogram type of prior distribution, [37]) can be obtained. This may complicate the posterior distribution in equation (10) and lead to difficulties in expressing it analytically and explicitly.As mentioned before,ξ¯ is a Gaussian random vector with a mean vector aμN+bl¯ and covariance matrix C¯. The likelihood function is thus given by(11)Pξ¯X,M=2π−n/2detC¯−1/2exp−12ξ¯−aμN+bl¯TC¯−1ξ¯−aμN+bl¯.The prior distribution and the likelihood function are then substituted into equation (10) to obtain the posterior distribution PXξ¯,M. The posterior distribution involves high dimension integration, which is one key limitation of Bayesian methods. To avoid the computational complexity, Markov Chain Monte Carlo Simulation (MCMCS) is adapted in this study to generate a large number of samples of random field parameters X from the posterior distribution (i.e., equation (10)). Note that MCMCS provides an appropriate tool to obtain the posterior knowledge in Bayesian analysis (e.g., [38]). Metropolis-Hastings (M-H) algorithm is used to generate a large number of X samples to describe the posterior distribution (i.e., equation (10), e.g., [15]). Samples of X are subsequently used to estimate the statistics of random field parameters through conventional statistical analyses.


## 4.2. Bayesian Selection of Correlation Function

The correlation functionM within a single statistically homogeneous clay layer is considered as a prescribed one with unknown λ in the previous section. This section presents a Bayesian model class selection method to determine the most probable correlation function M∗ among a pool of candidate correlation functions for a given set of ξ¯. Consider a number (i.e., NCF) of candidate correlation functions Mk, k = 1, 2, …, NCF within the clay layer, such as SECF, BNCF, SMCF, and SQECF shown in Figure 1. For a given set of ξ¯, the plausibility of the k-th candidate correlation function Mk in the clay layer is quantified by its occurrence probability PMkξ¯. Using Bayes’ theorem, PMkξ¯ is given by (e.g., [15, 36])(12)PMkξ¯=Pξ¯MkPMkPξ¯,k=1,2,…,NCF,where Pξ¯ is a normalizing constant that is independent of Mk and it is calculated as Pξ¯=∑k=1NCFPξ¯MkPMk; PMk is the prior probability of Mk; and Pξ¯Mk is the conditional probability density function of ξ¯ for a given correlation function Mk, and it quantifies the information on Mk provided by CPT data and is referred to as “evidence” for Mk provided by ξ¯.In the case of no prevailing prior knowledge, theNCF candidate correlation functions in the clay layer are considered to have the same prior probability PMk=1/NCF. Since Pξ¯ and PMk in equation (12) are constants, PMkξ¯ is proportional to Pξ¯Mk. The most probable correlation function M∗ is the one that has the maximum value of PMkξ¯, and it also has the maximum value of Pξ¯Mk. Note that Pξ¯Mk is actually the normalizing constant K in equation (10) by setting M as Mk, which involves a multidimensional integration. As discussed in the previous section, there is no need to calculate K (i.e., Pξ¯Mk) when MCMCS is adapted to the posterior distribution in equation (10) for Mk. However, it is needed herein for Bayesian selection of correlation function. MCMCS samples of X generated from equation (10) for Mk and Gaussian copula are used to calculate the evidence Pξ¯Mk (e.g., [15]), maintaining the advantage of no need to directly calculate the multidimensional integration (i.e., K).Using the Theorem of Total Probability, the evidencePξ¯Mk for Mk is expressed as(13)Pξ¯Mk=∫Pξ¯X,MkPXMkdX,where Pξ¯X,Mk is the likelihood function by setting M = Mk; and PXMk is the prior distribution of X. Then, the logarithm of evidence (i.e., lnPξ¯Mk) can be written as (e.g., [36])(14)lnPξ¯Mk=∫lnPξ¯X,MkPXξ¯,MkdX−∫lnPXξ¯,MkPXMkPXξ¯,MkdX.Based on theNp random samples of X estimated from PXξ¯,Mk, lnPξ¯Mk can be estimated as(15)lnPξ¯Mk≈∑i=1NplnPξ¯xi,MkNp−∑i=1NplnPxiξ¯,Mk/PxiMkNp,in which PxiMk is the prior distribution evaluated at i-th sample xi = [μi, σi, λi] of X, which is assumed to be a constant; Pξ¯xi,Mk is the likelihood function calculated using equation (11) for each xi; Pxiξ¯,Mk is posterior distribution evaluated using equation (10) at i-th sample xi.EvaluatingPxiξ¯,Mk needs the posterior distribution PXξ¯,Mk in equation (10), which is unknown and is numerically represented by Np MCMCS samples of X. To evaluate PXξ¯,Mk, copula theory is used to reconstruct PXξ¯,Mk based on Np samples of X in this paper. Details of calculation of the evidence using MCMCS and Gaussian copula can be found in Tian et al. [15].For a givenMk, the MCMCS-based Bayesian approach in the previous subsection is applied to generate random samples of X from the posterior distribution. These samples are subsequently used in equation (15) to evaluate lnPξ¯Mk. The procedure is repetitively performed NCF times for Mk, k = 1, 2, …, NCF to obtain their respective values of lnPξ¯Mk. Then, the most probable correlation function M∗ is determined by comparing the values of lnPξ¯Mk in each candidate correlation function. The one with the maximum value of lnPξ¯Mk is taken as M∗. Its corresponding posterior statistics (e.g., posterior mean values μ′, σ′, and λ′) estimated from MCMCS samples constitute the posterior knowledge on random field parameters. In other words, the most probable correlation function M∗ and posterior estimates (e.g., μ′, σ′, and λ′) of random field parameters are determined simultaneously for probabilistic characterization of the ISV of undrained shear strength based on both prior knowledge and CPT data.


## 5. Implementation Procedure

The Bayesian approach for probabilistic characterization comprises ten steps. Details of each step and its associated equations are summarized as follows:(1)
Obtain a set of cone tip resistanceqc and sleeve friction fs versus depth data for a clay layer from CPT, and convert them to the values of normalized cone resistance Qt and normalized friction ratio FR at different depths using equations (2) and (3), respectively.(2)
Calculate the soil behavior classification indexIc using equation (1) and identify homogeneous soil layers by examining whether the Cov of Ic is less than 0.1.(3)
Convert the values ofQt to ξ¯=ξD1,ξD2,…,ξDn at different depths D1, D2, …, Dn and perform conventional statistical analyses on the CPT data ξ¯ (i.e., ln Qt) to estimate its scale of fluctuation by fitting a candidate correlation function (e.g., SECF, BNCF, SMCF, and SQECF shown in Figure 1) to the sample autocorrelation function, calculate Bartlett statistic profiles using equation (4) where estimated scale of fluctuation lnQt is used, and compare the peak values of Bartlett statistic profiles Bmax with the critical values Bcrit to assess the stationarity of CPT data ξ¯.(4)
Choose the soil layers that satisfy both statistical homogeneity and weak stationarity to probabilistically characterize the ISV of undrained shear strength using random field model.(5)
DetermineNCF candidate correlation functions Mk, k = 1, 2, …, NCF, of a statistically homogeneous clay layer, such as SECF, BNCF, SMCF, and SQECF shown in Figure 1.(6)
Use M-H algorithm to generateNp samples of μ, σ, and λ based on equation (10) for a given Mk.(7)
Estimate the respective posterior statistics (e.g.,μ′, σ′, and λ′) based on their MCMCS samples through conventional statistical analyses.(8)
CalculatelnPξ¯Mk of Mk using equation (15), and the evidence Pξ¯Mk is then obtained by taking the exponential of lnPξ¯Mk.(9)
Repeat steps (6)–(8) NCF times to generate Np MCMCS samples of μ, σ, and λ and calculate lnPξ¯Mk and Pξ¯Mk for each candidate correlation function in the statistically homogeneous clay layer.(10)
Compare the values oflnPξ¯Mk (or Pξ¯Mk) of the NCF candidate correlation functions and select the one with the maximum value of lnPξ¯Mk as M∗ in the statistically homogeneous clay layer. Its corresponding posterior statistics and distributions reflect the posterior knowledge on μ, σ, and λ and probabilistically characterize the ISV of undrained shear strength of the statistically homogeneous clay layer.The ten steps can be readily implemented in commercial software packages, such as MATLAB (e.g., [39]). For geotechnical practitioners, the user function can be used as a “black box.” They need to provide project-specific test results (e.g., CPT data), prior knowledge (e.g., reasonable ranges of soil properties), transformation model, candidate correlation functions, and the number of MCMCS samples as input. The user function returns the most probable correlation function and its corresponding posterior estimates (e.g., mean values and standard deviations) of random field parameters within the statistically homogeneous layers at the specific site as output. This significantly improves the practicality of the Bayesian approaches. Identification of statistically homogeneous units and the Bayesian approach and its implementation procedure described above are illustrated through a set of real-life CPT data in the next section.


## 6. Illustrative Example

The Bayesian approaches are illustrated for probabilistic characterization of the ISV of undrained shear strength of soft clay at Yili station in Shanghai, China. Yili station serves as an interchange station for two orthogonal crossing subway lines in Shanghai (e.g., [40]). The site is comprised of six layers below the ground surface at a depth of about 49.2 m. The uppermost layer is a thin fill soil extended from the ground surface to about 0.8 m. The second layer, termed medium-soft clay, is slightly overconsolidated owing to desiccation. The third layer, the fourth layer, and the fifth layer are generally normally consolidated clay, with low permeability, high compressibility, and low strength. The sixth layer consists of medium compressible and slightly overconsolidated clay. The thicknesses of the six layers are 0.8 m, 2.7 m, 4.5 m, 9 m, 10.3 m, and 21.9 m, respectively. The groundwater table is located at about 1.5 m below the ground surface, and the total unit weight of the six layers is about 17.8 kN/m3. Figure 2 shows a set of cone tip resistance qc, sleeve friction fs, and pore water pressure u versus depth, respectively, obtained from CPT at this clay site. The average measurement interval is about 0.1 m. The profile of vertical effective stress, σv0′, is shown in Figure 2. The relationship between σv0′ and depth D is a linear equation: σv0′=8D+14.7 (kPa), where D is the depth below the ground surface varying from 1.5 m to 50 m. Figure 2 also shows the Su values indirectly obtained by equation (6) using CPT data and r values using the Su values and their corresponding σv0′. The CPT data is assessed for its statistical homogeneity and weak stationarity. Then, it is applied to identify random field model of r using the proposed Bayesian approaches.Figure 2
Soil profiles and CPT data at the Yili station.Consider, for example, four candidate correlation functions shown in Figure1, that is, SECF, BNCF, SMCF, and SQECF, which are denoted as M1, M2, M3, and M4 in this study, respectively. The prior knowledge on μ, σ, and λ is taken as a joint uniform distribution, which is equal to a constant (e.g., [36]). Based on the set of CPT data shown in Figure 2 and the prior knowledge, an MCMCS run is performed to generate 100,000 samples from the posterior distribution for each candidate correlation function. Note that the number of MCMCS samples is sufficient for MCMCS to reach the stationary condition of Markov Chain and to generate reasonably accurate probabilistic estimates of μ, σ, and λ in this study.


### 6.1. Results of Identification of Statistically Homogeneous Layers

Figure3 shows the normalized cone resistance Qt, normalized friction ratio FR, and soil behavior classification index Ic at Yili station. It can be seen that Ic ranges from 2.167 to 3.712, and the soil types cover a relatively wide range (zones 3 to 5), that is, from clay to sand mixtures. Figure 3 also shows that Ic of the medium clay layer and medium silty clay layer do not exhibit gradual obvious changes with depth. The mean values of Ic of the medium clay layer and medium silty clay layer are 0.02 and 0.03, respectively, which in no case exceed 0.1. Therefore, the medium clay layer and medium silty clay layer are essentially identified as homogeneous soil units and selected for subsequent assessment of weak stationarity and statistical analyses.Figure 3
Profiles of normalized cone resistanceQt, normalized friction ratio FR, and soil behavior classification index Ic at Yili station.
(a)
(b)
(c)The homogeneous soil units (i.e., medium clay layer and medium silty clay layer) at this site are further screened using the modified Bartlett test, which is capable of rejecting the null hypothesis of weak stationarity for spatially correlated data. The MBS procedure proposed by Phoon et al. [1] is used in this study to assess the weak stationarity of CPT data ξ¯. According to equation (7), there is a linear empirical relationship between ln r and CPT data ξ¯. If the profiles of ξ¯ have weak stationarity, r is also assumed as weak stationarity. This assumption would be validated using the Bayesian identified random field model of r. The scale of fluctuation of ξ¯ is estimated by fitting a candidate correlation function, as shown in Figure 1, to the initial part of the sample autocorrelation function. Then, the Bartlett statistic profiles of the homogeneous soil units (i.e., medium clay layer and medium silty clay layer) are calculated by the estimated scale of fluctuation.Figures4 and 5 show the Bartlett statistic profiles of the medium clay layer and medium silty clay layer. It is shown that the peak values of the Bartlett statistic profiles, Bmax, of these two clay layers are less than the critical values, Bcrit, respectively, which proves that the CPT data profiles of the medium clay layer and medium silty clay layer have weak stationarity. The CPT data ξ¯ obtained from these two clay layers can be used as input in the subsequent statistical analyses for probabilistic characterization of undrained shear strength. Then, under the Bayesian framework, CPT data and prior information are integrated to determine the random field parameters and select the most proper correlation function among a pool of candidates to characterize the ISV of r. It should be noted that the Bartlett statistic profiles of the medium clay layer and medium silty clay layer are calculated by equation (4), in which the scale of fluctuation of ξ¯ is robustly estimated by conventional statistical analyses.Figure 4
Modified Bartlett statistic profiles of CPT data for the medium clay layer. (a)M1. (b) M2. (c) M3. (d) M4.
(a)
(b)
(c)
(d)Figure 5
Modified Bartlett statistic profiles of CPT data for the medium silty clay layer. (a)M1. (b) M2. (c) M3. (d) M4.
(a)
(b)
(c)
(d)


### 6.2. Results of Random Field Model of Undrained Shear Strength Ratio

Table3 summarizes the results of random field parameters and correlation function of r of the medium clay layer and medium silty clay layer identified from the Bayesian approaches using the given set of CPT data (as shown in Figure 2) and prior knowledge. It is shown that, for the medium clay layer, M4 (i.e., SQECF) has the largest value of PMkξ¯ (i.e., 0.293) among M1–M4. Therefore, it is taken as M∗ for probabilistically characterizing r of the medium clay layer at Yili station. Similar observation is also obtained for the medium silty clay layer. M4 (i.e., SQECF) is also taken as M∗ for probabilistically characterizing r of the medium silty clay layer. Table 3 also includes the respective posterior estimates (e.g., μ′, σ′, and λ′) of μ, σ, and λ of the medium clay layer and medium silty clay layer for M1–M4. The posterior estimates of μ corresponding to M4 (i.e., M∗) of the medium clay layer and medium silty clay layer are 0.21 and 0.24, respectively. These values fall within the reported values (e.g., [10, 16, 41]. The posterior Cov of r corresponding to M4 (i.e., M∗) of the two clay layers are 9.5% and 16.7%, respectively. The estimated values of Cov also essentially fall within the reported values (i.e., 10% to 35%).Table 3
Results of random field parameters and correlation function ofr of the medium clay layer and medium silty clay layer obtained from the Bayesian approaches.
Soil type
Correlation function
P
M
k
ξ
¯
μ
′
σ
′
λ
′ (m)
Medium clay layer
M
1
0.251
0.21
0.01
3.41
M
2
0.189
0.21
0.01
3.33
M
3
0.267
0.21
0.01
3.58
M
4
a
0.293
0.21
0.02
3.44
Medium silty clay layer
M
1
0.105
0.23
0.03
4.41
M
2
0.173
0.23
0.03
4.04
M
3
0.281
0.23
0.03
4.46
M
4
a
0.441
0.24
0.04
4.30
aMost probable correlation function.


### 6.3. Comparison with Conventional Statistical Analyses

For comparison, conventional statistical analyses are performed on ther profiles (as shown in Figure 2) estimated from CPT data using the empirical regression function shown in equation (6). μ, σ, and sample autocorrelation function are directly estimated from the r profiles. λ and M∗ are obtained by fitting a candidate correlation function (e.g., M1–M4) to the sample autocorrelation function and comparing the goodness-of-fitting in different candidate correlation functions (e.g., [20, 21]). The results of estimates of μ, σ, λ, and M∗ are compared with those obtained from the Bayesian approaches.Table4 summarizes the results estimated from the conventional statistical analyses, including the direct estimates of random field parameters and goodness-of-fitting for each candidate correlation function (e.g., M1–M4). Table 4 shows that M1 (i.e., SECF) is selected as M∗ among M1–M4 for the medium clay layer and medium silty clay layer, which is obviously different from M∗ determined by the Bayesian approaches (see Table 3). Table 4 also includes the direct estimates of random field parameters (see columns 4-5 in Table 4) of the medium clay layer and medium silty clay layer, respectively. Compared with the results obtained from the Bayesian approaches (see Table 3), the differences between μ (i.e., 0.01) estimated from these two approaches are relatively small. However, the σ values (see column 5 in Table 4) estimated from conventional statistical analyses are greater than those from the Bayesian approaches (see column 5 in Table 3), and the λ values (see column 6 in Table 4) obtained from conventional statistical analyses are less than the Bayesian approaches (see column 6 in Table 3). Such a difference might be attributed to the reason that the effect of transformation uncertainty is not explicitly considered in conventional statistical analyses (e.g., [15]). Both the ISV of γ and the transformation uncertainty propagate into the profiles of r directly estimated from the CPT data using the regression function (i.e., equation (6)). Therefore, the sample autocorrelation function estimated from the r profiles does not reflect the actual spatial correlation. The variability in the r profiles contains not only the ISV but also the transformation uncertainty, and it is, therefore, greater than the actual ISV of r, which subsequently leads to overestimation in σ and underestimation in λ. It is also not surprising to see that the correlation function most directly obtained by fitting the sample autocorrelation function of the r profiles is inconsistent with that estimated from the Bayesian approaches. In contrast, the transformation uncertainty is, explicitly and rationally, taken into account in the Bayesian approaches.Table 4
Results of random field parameters and correlation function ofr of the medium clay layer and medium silty clay layer obtained from conventional statistical analyses.
Soil type
Correlation function
Goodness-of-fittingR2
μ
′
σ
′
λ
′ (m)
Medium clay layer
M
1
∗
0.848
0.21
0.02
1.27
M
2
0.722
0.21
0.02
1.34
M
3
0.767
0.21
0.02
1.15
M
4
0.722
0.21
0.02
1.08
Medium silty clay layer
M
1
∗
0.906
0.23
0.05
0.79
M
2
0.822
0.23
0.05
0.63
M
3
0.844
0.23
0.05
0.67
M
4
0.810
0.23
0.05
0.59
∗Best-fit correlation function.


### 6.4. Validation of Bayesian Identification of Random Field Model

For validation, the most probable correlation functions ofr of the medium clay layer and medium silty clay layer obtained from the Bayesian approaches are compared with the sample autocorrelation functions, as shown in Figures 6(a) and 6(b). It is shown that the initial parts of circles plot close to the solid lines, indicating that the proposed Bayesian approaches properly identify random field model for probabilistic characterization of ISV of undrained shear strength ratio r.Figure 6
Comparison of correlation functions estimated from the Bayesian approaches and the sample autocorrelation function. (a) Medium clay layer. (b) Medium silty clay layer.
(a)
(b)It should be noted that the weak stationarity of CPT dataξ¯ is assessed by the scale of fluctuation estimated by conventional statistical analyses (i.e., estimated by fitting a candidate correlation function, as shown in Figure 1, to the sample autocorrelation function). For illustration, the scale of fluctuation of ξ¯ obtained from the Bayesian approaches is used to validate the weak stationarity of CPT data within the medium clay layer and medium silty clay layer. Figures 7 and 8 show the Bartlett statistic profiles of CPT data ξ¯ in the medium clay layer and medium silty clay layer are calculated by the scale of fluctuation obtained from the Bayesian approaches. It can be seen that, for the medium clay layer and medium silty clay layer, the peak values of the Bartlett statistic profiles, Bmax, are less than the critical values, Bcrit. In Section 6.1, the r profiles of the medium clay layer and medium silty clay layer are assumed to have stationarity if the CPT data profiles have weak stationarity. This assumption is further verified using the random field parameters of r obtained from Bayesian approaches (as shown in Table 3). Figures 9 and 10 show that the peak values Bmax of r are less than the critical values, Bcrit. It is suggested that the hypothesis of stationarity of r profiles is accepted and the random field parameters and correlation function of r obtained from the Bayesian approaches can be used to properly characterize the ISV of r.Figure 7
Validation of stationarity of CPT data for the medium clay layer. (a)M1. (b) M2. (c) M3. (d) M4.
(a)
(b)
(c)
(d)Figure 8
Validation of stationarity of CPT data for the medium silty clay layer. (a)M1. (b) M2. (c) M3. (d) M4.
(a)
(b)
(c)
(d)Figure 9
Validation of stationarity ofr for the medium clay layer. (a) M1. (b) M2. (c) M3. (d) M4.
(a)
(b)
(c)
(d)Figure 10
Validation of stationarity ofr for the medium silty clay layer. (a) M1. (b) M2. (c) M3. (d) M4.
(a)
(b)
(c)
(d)


## 6.1. Results of Identification of Statistically Homogeneous Layers

Figure3 shows the normalized cone resistance Qt, normalized friction ratio FR, and soil behavior classification index Ic at Yili station. It can be seen that Ic ranges from 2.167 to 3.712, and the soil types cover a relatively wide range (zones 3 to 5), that is, from clay to sand mixtures. Figure 3 also shows that Ic of the medium clay layer and medium silty clay layer do not exhibit gradual obvious changes with depth. The mean values of Ic of the medium clay layer and medium silty clay layer are 0.02 and 0.03, respectively, which in no case exceed 0.1. Therefore, the medium clay layer and medium silty clay layer are essentially identified as homogeneous soil units and selected for subsequent assessment of weak stationarity and statistical analyses.Figure 3
Profiles of normalized cone resistanceQt, normalized friction ratio FR, and soil behavior classification index Ic at Yili station.
(a)
(b)
(c)The homogeneous soil units (i.e., medium clay layer and medium silty clay layer) at this site are further screened using the modified Bartlett test, which is capable of rejecting the null hypothesis of weak stationarity for spatially correlated data. The MBS procedure proposed by Phoon et al. [1] is used in this study to assess the weak stationarity of CPT data ξ¯. According to equation (7), there is a linear empirical relationship between ln r and CPT data ξ¯. If the profiles of ξ¯ have weak stationarity, r is also assumed as weak stationarity. This assumption would be validated using the Bayesian identified random field model of r. The scale of fluctuation of ξ¯ is estimated by fitting a candidate correlation function, as shown in Figure 1, to the initial part of the sample autocorrelation function. Then, the Bartlett statistic profiles of the homogeneous soil units (i.e., medium clay layer and medium silty clay layer) are calculated by the estimated scale of fluctuation.Figures4 and 5 show the Bartlett statistic profiles of the medium clay layer and medium silty clay layer. It is shown that the peak values of the Bartlett statistic profiles, Bmax, of these two clay layers are less than the critical values, Bcrit, respectively, which proves that the CPT data profiles of the medium clay layer and medium silty clay layer have weak stationarity. The CPT data ξ¯ obtained from these two clay layers can be used as input in the subsequent statistical analyses for probabilistic characterization of undrained shear strength. Then, under the Bayesian framework, CPT data and prior information are integrated to determine the random field parameters and select the most proper correlation function among a pool of candidates to characterize the ISV of r. It should be noted that the Bartlett statistic profiles of the medium clay layer and medium silty clay layer are calculated by equation (4), in which the scale of fluctuation of ξ¯ is robustly estimated by conventional statistical analyses.Figure 4
Modified Bartlett statistic profiles of CPT data for the medium clay layer. (a)M1. (b) M2. (c) M3. (d) M4.
(a)
(b)
(c)
(d)Figure 5
Modified Bartlett statistic profiles of CPT data for the medium silty clay layer. (a)M1. (b) M2. (c) M3. (d) M4.
(a)
(b)
(c)
(d)


## 6.2. Results of Random Field Model of Undrained Shear Strength Ratio

Table3 summarizes the results of random field parameters and correlation function of r of the medium clay layer and medium silty clay layer identified from the Bayesian approaches using the given set of CPT data (as shown in Figure 2) and prior knowledge. It is shown that, for the medium clay layer, M4 (i.e., SQECF) has the largest value of PMkξ¯ (i.e., 0.293) among M1–M4. Therefore, it is taken as M∗ for probabilistically characterizing r of the medium clay layer at Yili station. Similar observation is also obtained for the medium silty clay layer. M4 (i.e., SQECF) is also taken as M∗ for probabilistically characterizing r of the medium silty clay layer. Table 3 also includes the respective posterior estimates (e.g., μ′, σ′, and λ′) of μ, σ, and λ of the medium clay layer and medium silty clay layer for M1–M4. The posterior estimates of μ corresponding to M4 (i.e., M∗) of the medium clay layer and medium silty clay layer are 0.21 and 0.24, respectively. These values fall within the reported values (e.g., [10, 16, 41]. The posterior Cov of r corresponding to M4 (i.e., M∗) of the two clay layers are 9.5% and 16.7%, respectively. The estimated values of Cov also essentially fall within the reported values (i.e., 10% to 35%).Table 3
Results of random field parameters and correlation function ofr of the medium clay layer and medium silty clay layer obtained from the Bayesian approaches.
Soil type
Correlation function
P
M
k
ξ
¯
μ
′
σ
′
λ
′ (m)
Medium clay layer
M
1
0.251
0.21
0.01
3.41
M
2
0.189
0.21
0.01
3.33
M
3
0.267
0.21
0.01
3.58
M
4
a
0.293
0.21
0.02
3.44
Medium silty clay layer
M
1
0.105
0.23
0.03
4.41
M
2
0.173
0.23
0.03
4.04
M
3
0.281
0.23
0.03
4.46
M
4
a
0.441
0.24
0.04
4.30
aMost probable correlation function.


## 6.3. Comparison with Conventional Statistical Analyses

For comparison, conventional statistical analyses are performed on ther profiles (as shown in Figure 2) estimated from CPT data using the empirical regression function shown in equation (6). μ, σ, and sample autocorrelation function are directly estimated from the r profiles. λ and M∗ are obtained by fitting a candidate correlation function (e.g., M1–M4) to the sample autocorrelation function and comparing the goodness-of-fitting in different candidate correlation functions (e.g., [20, 21]). The results of estimates of μ, σ, λ, and M∗ are compared with those obtained from the Bayesian approaches.Table4 summarizes the results estimated from the conventional statistical analyses, including the direct estimates of random field parameters and goodness-of-fitting for each candidate correlation function (e.g., M1–M4). Table 4 shows that M1 (i.e., SECF) is selected as M∗ among M1–M4 for the medium clay layer and medium silty clay layer, which is obviously different from M∗ determined by the Bayesian approaches (see Table 3). Table 4 also includes the direct estimates of random field parameters (see columns 4-5 in Table 4) of the medium clay layer and medium silty clay layer, respectively. Compared with the results obtained from the Bayesian approaches (see Table 3), the differences between μ (i.e., 0.01) estimated from these two approaches are relatively small. However, the σ values (see column 5 in Table 4) estimated from conventional statistical analyses are greater than those from the Bayesian approaches (see column 5 in Table 3), and the λ values (see column 6 in Table 4) obtained from conventional statistical analyses are less than the Bayesian approaches (see column 6 in Table 3). Such a difference might be attributed to the reason that the effect of transformation uncertainty is not explicitly considered in conventional statistical analyses (e.g., [15]). Both the ISV of γ and the transformation uncertainty propagate into the profiles of r directly estimated from the CPT data using the regression function (i.e., equation (6)). Therefore, the sample autocorrelation function estimated from the r profiles does not reflect the actual spatial correlation. The variability in the r profiles contains not only the ISV but also the transformation uncertainty, and it is, therefore, greater than the actual ISV of r, which subsequently leads to overestimation in σ and underestimation in λ. It is also not surprising to see that the correlation function most directly obtained by fitting the sample autocorrelation function of the r profiles is inconsistent with that estimated from the Bayesian approaches. In contrast, the transformation uncertainty is, explicitly and rationally, taken into account in the Bayesian approaches.Table 4
Results of random field parameters and correlation function ofr of the medium clay layer and medium silty clay layer obtained from conventional statistical analyses.
Soil type
Correlation function
Goodness-of-fittingR2
μ
′
σ
′
λ
′ (m)
Medium clay layer
M
1
∗
0.848
0.21
0.02
1.27
M
2
0.722
0.21
0.02
1.34
M
3
0.767
0.21
0.02
1.15
M
4
0.722
0.21
0.02
1.08
Medium silty clay layer
M
1
∗
0.906
0.23
0.05
0.79
M
2
0.822
0.23
0.05
0.63
M
3
0.844
0.23
0.05
0.67
M
4
0.810
0.23
0.05
0.59
∗Best-fit correlation function.


## 6.4. Validation of Bayesian Identification of Random Field Model

For validation, the most probable correlation functions ofr of the medium clay layer and medium silty clay layer obtained from the Bayesian approaches are compared with the sample autocorrelation functions, as shown in Figures 6(a) and 6(b). It is shown that the initial parts of circles plot close to the solid lines, indicating that the proposed Bayesian approaches properly identify random field model for probabilistic characterization of ISV of undrained shear strength ratio r.Figure 6
Comparison of correlation functions estimated from the Bayesian approaches and the sample autocorrelation function. (a) Medium clay layer. (b) Medium silty clay layer.
(a)
(b)It should be noted that the weak stationarity of CPT dataξ¯ is assessed by the scale of fluctuation estimated by conventional statistical analyses (i.e., estimated by fitting a candidate correlation function, as shown in Figure 1, to the sample autocorrelation function). For illustration, the scale of fluctuation of ξ¯ obtained from the Bayesian approaches is used to validate the weak stationarity of CPT data within the medium clay layer and medium silty clay layer. Figures 7 and 8 show the Bartlett statistic profiles of CPT data ξ¯ in the medium clay layer and medium silty clay layer are calculated by the scale of fluctuation obtained from the Bayesian approaches. It can be seen that, for the medium clay layer and medium silty clay layer, the peak values of the Bartlett statistic profiles, Bmax, are less than the critical values, Bcrit. In Section 6.1, the r profiles of the medium clay layer and medium silty clay layer are assumed to have stationarity if the CPT data profiles have weak stationarity. This assumption is further verified using the random field parameters of r obtained from Bayesian approaches (as shown in Table 3). Figures 9 and 10 show that the peak values Bmax of r are less than the critical values, Bcrit. It is suggested that the hypothesis of stationarity of r profiles is accepted and the random field parameters and correlation function of r obtained from the Bayesian approaches can be used to properly characterize the ISV of r.Figure 7
Validation of stationarity of CPT data for the medium clay layer. (a)M1. (b) M2. (c) M3. (d) M4.
(a)
(b)
(c)
(d)Figure 8
Validation of stationarity of CPT data for the medium silty clay layer. (a)M1. (b) M2. (c) M3. (d) M4.
(a)
(b)
(c)
(d)Figure 9
Validation of stationarity ofr for the medium clay layer. (a) M1. (b) M2. (c) M3. (d) M4.
(a)
(b)
(c)
(d)Figure 10
Validation of stationarity ofr for the medium silty clay layer. (a) M1. (b) M2. (c) M3. (d) M4.
(a)
(b)
(c)
(d)


## 7. Summary and Conclusions

This paper presented Bayesian approaches for probabilistic characterization of ISV of undrained shear strength within a statistically homogeneous clay layer based on CPT data and prior knowledge. Homogeneous soil units were first identified from the CPT data using the soil behavior classification index. The identified homogeneous soil units were further assessed using modified Bartlett’s test to reject the null hypothesis of weak stationarity. Then, the Bayesian approaches were used to identify random field parameters (i.e., meanμ, standard deviation σ, and scale of fluctuation λ) and simultaneously select the most probable correlation function among a pool of candidate correlation functions using identified homogeneous or stationary CPT profiles and prior knowledge. MCMCS and Gaussian copula-based method were used to bypass multidimensional integration involved in the Bayesian approaches.Equations were derived for identification of statistically homogeneous soil units and the Bayesian approaches. They were illustrated using profiles of CPT data obtained from a clay site in Shanghai, China. The statistically homogeneous clay layers were identified from the site-specific CPT data. Random field parameters and correlation function in a statistically homogeneous clay layer obtained from the Bayesian approaches were compared with those estimated from conventional statistical analyses. It has been shown that the variability in the profiles of undrained shear strength directly estimated from the CPT data using the empirical regression function contains not only the ISV of undrained shear strength but also the transformation uncertainty, which subsequently leads to overestimation inσ and underestimation in λ. The Bayesian approaches provide a rational tool to properly identify random field parameters and correlation function for probabilistic characterization of ISV of undrained shear strength using CPT data with proper consideration of transformation uncertainty.In addition, the results were validated by comparing the sample autocorrelation function with correlation function obtained from the Bayesian approaches. The assessment of stationarity was also validated based on the random field parameters obtained from Bayesian approaches. It was shown that the most probable correlation function obtained from the Bayesian approaches agrees well with the sample autocorrelation function. The hypothesis of stationarity of CPT data and undrained shear strength ratio profiles is accepted. The random field parameters and correlation function of undrained shear strength ratio,r, obtained from the Bayesian approaches can be used to properly characterize the ISV of r.



---

*Source: 9617698-2020-02-07.xml*