Uncertainty Analysis of Creep Behavior of Compacted Loess and a Non-Deterministic Predication Method for Post-Construction Settlement of a High-Fill Embankment

Abstract

Property of geotechnical materials has inherent uncertainty due to the complex formation process and inevitable test error. However, existing long-term deformation prediction methods for geotechnical structure such as a filling embankment are deterministic, which ignores the uncertainty of soil property. In this study, the uncertainty of creep behavior of compacted loess was investigated through repetitive creep tests and statistical analysis. Five different loading levels and two loading modes were considered in the tests. The creep test was repeated 45 times for each condition. Through a statistical analysis for the test results, a modified Merchant creep model was established to improve the accuracy of long-term deformation prediction. An empirical transformation equation between staged loading and separated loading mode of the creep test results was also introduced to improve applicability of the method. On this basis, a non-deterministic predication method for post-construction settlement of loess fill embankment was proposed. Furthermore, the proposed method was applied to the prediction of the post-construction of a 61.5 m loess filling embankment. The measured on-site post-construction settlement value falls within the 95% confidence interval of the predicted range which proves the efficiency and practicability of the proposed non-deterministic predication method. Compared to deterministic methods, the proposed method can describe the predicted deformation in a probabilistic way in the form of contour plot. The proposed method provides a basic approach for the probabilistic design and reliability assessment of filling engineering.

1. Introduction

Accurate deformation prediction is the premise of safety and quality of geotechnical engineering [1]. Post-construction settlement, a macro-outcome of the creep behavior of filled soil, is the main influencing factor of subsequent construction of geo-structures such as high fill embankments [2]. To predict the post-construction accurately so as to guide the subsequent construction and design locating on the filling region [3], creep behavior of the filled soil must be studied thoroughly. Creep tests are the data source for the soil creep model, which has two major forms, the consolidation creep tests [4] and the triaxial creep tests [5]. The former one is for a lateral limiting condition such as settlement of soil foundation, while the latter one is for a triaxial stress state such as the long-term slip deformation of a filling slope. Based on creep tests, a soil creep model can be formed. For example, Yin et al. [6] proposed a fractal derivative viscoelastic plastic creep model considering the effect of damage which requires fewer parameters and has higher simulation precision. Ausilio and Conte [7] introduced a simplified equation to link the settlement rate of unsaturated soil with the average degree of consolidation and analyzed the one-dimensional consolidation of unsaturated soil. Conte [8] carried out consolidation analysis of unsaturated soil for plane strain and axisymmetric problems. They also extended the analysis results to the consolidation of coupled and uncoupled unsaturated soil. The effects of water, air coupling and uncoupling on the consolidation characteristics of unsaturated soil were also compared and analyzed. Much in-depth research regarding creep models was also carried out [9,10,11,12,13,14,15,16]. To predict the post construction settlement, creep models were often applied to a numerical method such as the finite element method (FEM) [17], artificial neural network (ANN) [18], discrete element method (DEM) [19], etc. Using different creep description models and analysis methods, long-term deformation of different engineering cases was investigated. Meng et al. [20] investigated the effect of preloading on post-construction consolidation settlement of soft clay subjected to repeated loading after removal of a part of the preload. Chung [19] proposed the FEM-DEM contact coupling algorithm that provides an effective way to determine the stress distribution in contacting structures and can be applied to solving solid–structure interaction problems. Zheng [21] improved the contact coupling algorithm, which was further verified by an example. The extrapolation method is another way of post-construction settlement prediction which uses acquired monitoring data to predict the subsequent settlement. Wan and Doherty [22] established a data-driven approach for forecasting the behavior of embankments on soft soil during and shortly after construction. The method has high precision but cannot be used when no monitoring data are available. This research investigates the prediction method of long-term deformation with different tools and methods, but they share one common defect that the creep behavior of soil was treated as constant, and the long-term deformation was calculated in a deterministic way, which ignores the natural randomness of soil.

Different from artificial materials such as fiber-reinforced polymer [23,24], soil has inherent uncertainty [25] which affects the deformation property. Randomness of the properties of geotechnical materials consists of two parts, its inherent randomness due to the deposition history or formation process [26] and the test errors influenced by different devices, test operators and sample preparation quality [27]. All these components are expressed in the form of statistical variability of the resulting parameters from repetitive laboratory tests in practice. Considering the randomness of soils, strength properties were studied, and reliability of geo-structures was investigated in much of the literature [28,29,30,31,32]. For example, Yuan et al. [33] studied the “curve local average” of a two-dimensional random field and proposed a simple reliability analysis method for homogeneous natural slope. Halder et al. [34] combined the finite difference method with the random field model to study the influence of the spatial variability and randomness of soil strength parameters on the load-settlement response of strip foundation on the slope. Regarding the deformation problems, Gong et al. [35] proposed a framework for the probabilistic analysis of tunnel longitudinal performance that considered a conditional random field, by which tunnel longitudinal performance was analyzed in a probabilistic way. Yu et al. [36] investigate tunnel liner performance based on the concept of reliability-based design. A probabilistic code that evaluates the tunnel liner performance is described. Puglia [37] combined the random field theory and lattice discrete element method and considered the randomness of material properties in the contact constitutive model in their research. As a result, a random discrete element method was realized. Mahmoud and Soheil [25] performed a probabilistic analysis of shallow foundation settlement using a probabilistic Monte Carlo simulation. These studies focus on the influence of probabilistic features of soil on reliability problems or short-term instantaneous deformation. On the other hand, influence of the random features of soils on the long-term deformation such as the post-construction settlement of filling embankment has not been investigated yet. It is of great significance to introduce randomness in a creep model for soil and establish a method that predicts the post-construction settlement in a non-deterministic way, so as to provide more comprehensive information for the construction and design of geo-structures such as filling embankment, which is the objective of this paper.

In this study, through repetitive consolidation creep tests, the creep deformation of compacted loess under staged and separate loading modes was tested and discussed. The Merchant creep model is modified to increase its accuracy in long-term deformation prediction for compacted loess. A transformation formula between the test results of staged and separate loading mode is deduced to increase the practicability of uncertainty analysis. On these bases, a non-deterministic prediction model considering randomness of the creep behavior of compacted loess is proposed based on the statistical parameter analysis. An engineering application was finally carried out to prove its validity.

The remainder of the paper is organized as follows. Section 2 describes the repetitive tests of the creep behavior of loess. Section 3 discussed the test results. Section 4 introduced a conversion method that unifies the test data of different loading modes of creep tests. In Section 5, a randomness analysis for creep model parameters based on the test data was conducted. In Section 6, a non-deterministic prediction method for post-construction settlement of loess high-fill embankments was proposed and verified through engineering application. Finally, in Section 7, conclusions are provided.

2. Materials and Test Method

The purpose of a non-deterministic analysis is to quantify the effect of uncertainty of input on the system outputs [38]. In a non-deterministic analysis, quantitative description of input uncertainty is the first and the most crucial step. To capture the statistical feature of the creep behavior of compacted loess, a group of repetitive creep tests was conducted.

2.1. Materials

The test soil sample is collected from the high fill project in the new campus of Yan’an University, Yan’an City, Shaanxi Province. The soil sample is composed of silt, containing a small amount of silty clay, which belongs to Q₃ loess. Basic geotechnical tests were conducted to obtain the basic properties of the soil samples.

Table 1. Indicators of basic physical properties of soil samples.

Specific Weight	Liquid Limit/%	Plastic Limit/%	Plasticity Index/%	Particle Composition/%
				>0.075 mm	0.075~0.005 mm	<0.005 mm
2.70	29.70	18.40	11.30	1.05%	78.43%	20.52%

shows the resulting physical and mechanical properties. Particle-Size Distribution analysis was conducted using the hydrometer method following test standard GB/T 50123 (published by Ministry of Housing and Urban Rural Development of China). The resulting grading curve is shown in Figure 1. In the tests of the present study, samples with different water content were prepared by water film transfer method for humidification and placed in a moisturizing cylinder for 48 h before being used [39].

2.2. Test Method

To study the statistical feature of the creep behavior of compacted loess and establish the stochastic creep model for non-deterministic predication, this paper designed and performed a set of repetitive creep tests. To increase the practicability of our study, a high pressure oedometer, which is widely used in engineering practice, was used to conduct the creep tests in this paper.

Staged loading mode and separated loading mode were both chosen to analyze the statistical difference between the loading modes. To ensure good contact between the sample and the loading system, 25 kPa of preloading was exerted until the deformation is stable before the formal test. Then, five load stages (100 kPa, 200 kPa, 400 kPa, 800 kPa and 1600 kPa) were applied to test samples. The stability standard of the test is that the vertical cumulative deformation is less than 0.005 mm within 24 h. To prevent water loss of the sample during the test, the surface of the sample was wrapped with plastic wrap [40]. Additionally, water content was rechecked after the creep test. Dry density and water content of the test samples were set to 1.68 g/cm³ and 10%, respectively, which are the maximum dry density and optimum water content according to the result of compaction test. The creep test was repeated forty-five times for statistical analysis. Numbering and settings of the creep tests were shown in

Table 2. Numbering and settings of creep tests.

Sample NO.	Loading Type	Loading Level	Repetitive Times	Dry Density/ (g/cm3)	Water Content (%)
C1-1~C1-45; C2-1~C2-45; C3-1~C3-45; C4-1~C4-45; C5-1~C5-45	Separated loading	100 kPa, 200 kPa, 400 kPa, 800 kPa, 1600 kPa	45	1.68	10%
S1-1~S1-45; S2-1~S2-45; S3-1~S3-45; S4-1~S4-45; S5-1~S5-45	Staged loading

. Photos of the test process are shown in Figure 2.

3. Test Results

Test results were statistically analyzed in this section. Furthermore, a modified merchant model will be established in Section 4, and a parameters randomness analysis will also be conducted to quantitatively descript the statistical feature of the creep behavior of compacted loess. Figure 3a,b shows the typical relationship curve between creep strain and time under two different loading modes. Figure 2a is the unprocessed-graded loading results, Figure 2b is the comparison of graded loading and separate loading which shows the difference between two loading mode. The conversion equation was proposed accordingly. It can be observed from the figures that under separated loading, creep deformation increases rapidly with the applied load and converges over time. Under staged loading, the creep deformation changes step by step with the applied loads. The samples all went through three stages: a transient deformation stage, deceleration creep stage and stable creep stage. The soil experienced structural damage and reorganization during the creep process. The initial structure was damaged due to the load, and creep deformation became stable due to the completion of soil reorganization at a later creep period. The Boltzmann superposition principle [41] was used to translate the results of the staged loading to a separated form for comparison as shown in Figure 3b. It can be observed that, at all loading levels except 100 kPa, creep deformation obtained by the separated loading mode is greater than those of the staged loading mode, especially at 1600 kPa. The error of the translation accumulated with the growth of the loading level. The translated curve tends to be bigger than the original result of separated loading because the staged loading takes more time and tends to produce bigger deformation as the loading level increases. Theoretically, the results of separated loading and translated staged loading at 100 kPa should be identical since no translation is performed at this level. The obtained separated loading result is slightly smaller which is merely caused by test error. Since both loading modes are commonly used in practical engineering, it is necessary to set up a bridge between the two modes so that diverse types of data can both be used in the statistical analysis. This work will be discussed in Section 4.

As for the repetitive results, Figure 3c shows that under different loading levels, creep deformation exhibited evident randomness. Loading level has a positive correlation with the randomness of creep deformation. This is not surprising since a greater magnitude tends to cause a bigger possible error. In other words, with the same error rate, data with bigger magnitude show more obvious fluctuation. It should be noted that the repetitive tests were conducted by the same skilled test operator under the guidance of the standard test procedure and on a well-functioning test device. However, the test results still differ greatly among different test samples. The test errors that cause this randomness are due to the nonuniformity of the disturbed soil, sample preparation error (including the error of water content and dry density), measurement error and manual error (no matter how skillful the tester is). These errors can be reduced by following the standard operating procedure or by using a more sophisticated device but can never be completely eliminated. Therefore, when using the test data to predict the possible creep deformation for practical engineering, this randomness must be considered, and a non-deterministic predication should be conducted instead a deterministic one. Moreover, repetitive intensity also influences the result, which should be considered in engineering application.

4. Conversion Correction between Staged Loading and Separated Loading

Separated loading and staged loading both have their own advantages and limitations. The advantage of separated loading is its accuracy of different loading levels since they were treated separately. However, different samples may cause bigger errors during the sample preparation and load application. The test complexity is also higher. On the other hand, the advantage of staged loading is that different loading levels were exerted on one sample so that the sample preparation error can be excluded. However, the translated curve may be distorted. These two loading modes were both commonly used in practical engineering; therefore, according to the average results of the repetitive tests using these two loading modes, a conversion correction method was proposed so that the results of both modes can be used in the statistical analysis.

Firstly, an empirical fitting model was employed to analyze the test results. The fitting process was performed with software ORIGIN (version 2019b). The least square method and Levenberg–Marquardt algorithm were used for the linear and power function fitting described below, respectively. Using the average data of the two loading modes, the relationship between $t/\left[\epsilon \right(t){]}^{\ast }$ and t was fitted using the linear function (Equation (1)). Fitting results are shown in Figure 4. To increase the fitting precision, extra tests considering the loading level at 300 kPa, 600 kPa and 1200 kPa were conducted and analyzed in this section. It can be seen from Figure 4 that the relationship between $t/\left[\epsilon \right(t){]}^{\ast }$ and t is clearly linear for all loading levels and both loading modes.

$$\frac{t}{[\epsilon (t){]}^{\ast }}=kt+b$$

(1)

where $[\epsilon \left(t\right){]}^{\ast }$ is the average creep strain of the repetitive test; k and b are the slope and intercept for the linear relation between $t/\left[\epsilon \right(t){]}^{\ast }$ and t.

Then, the empirical fitting model for the studied loess can be obtained by translating Equation (1) to Equation (2):

$$[\epsilon (t){]}^{\ast }=\frac{t}{kt+b}$$

(2)

According to Equation (2), creep strain at any time can be calculated. When $t\to \infty$ the ultimate creep strain can be obtained as:

$$\epsilon \ast (t\to \infty )=\frac{1}{k}$$

(3)

The relationship between the fitting parameters k and b with the loading P is shown in Figure 5a,b, respectively. Fitting results show a strong exponential relation (with correlation coefficient bigger than 0.985) between k, b and loading level P, indicating a convergence trend for the slope and intercept for $t/\left[\epsilon \right(t){]}^{\ast }$-t. The slope tends to converge at 25 when P is big enough (Figure 5a), while the intercept approaches 0 with the growth of P. Their relationship can be expressed as a power function as shown in Equation (4):

$$k=A{p}^B,b=C{p}^D$$

(4)

where A, B, C, D are fitting parameters. Substituting Equation (4) for Equation (2), the empirical creep model of the studied loess can be obtained as:

$$\epsilon (t)=\frac{t}{A{p}^{B}t+C{p}^D}$$

(5)

Based on the obtained empirical creep model, slopes of the relation curve of $t/\left[\epsilon \right(t){]}^{\ast }$ and t for separated loading and staged loading were denoted by A_b and A_j. The relationship between A_b and A_j under different loading levels was shown in Figure 6. The relationship between A_b and A_j can be expressed as a power function as follows:

$$A_b{=f(A}_j{), A}_b=0.4278A_j{}^{1.1744}$$

(6)

According to Equation (3), when $t\to \infty$:

$${\epsilon }_b\ast (t\to \infty )=\frac{1}{A_b}{,\epsilon }_j\ast (t\to \infty )=\frac{1}{A_j}$$

(7)

Combining Equations (6) and (7), and generalizing the translation of the ultimate value to the whole curve, the following translation equation from staged loading to separated loading can be obtained:

$${\epsilon }_b=\frac{1}{f(\frac{1}{{\epsilon }_j})}{, \epsilon }_b=\frac{1}{0.4278}{\epsilon }_j{}^{1.1774}$$

(8)

where ${\epsilon }_j$ is the creep strain of staged loading mode while ${\epsilon }_b$ is that of separated loading mode. Using Equation (8), the creep test results from staged loading can be translated into separated loading.

5. Randomness Analysis for Creep Model Parameters

The original Merchant model is a three-parameter component model formed by the Hooke model and Kelvin model in series. The equation of original Merchant model is shown in Equation (9):

$$\epsilon (t)=\frac{{\sigma }_0}{E_0}+\frac{{\sigma }_0}{E_1}\left[1-e^{\left(-\frac{E_1}{{\eta }_1}t\right)}\right]$$

(9)

where $\epsilon \left(t\right)$ is the creep strain at time t. ${\sigma }_0$ is the vertical loading. $E_0$ is the coefficient of elasticity for Hooke model. $E_1$ and ${\eta }_1$ are the coefficient of elasticity and viscosity for Kelvin model, respectively.

The Merchant model cannot precisely describe the creep deformation from the decay creep stage to the stable creep stage. To address this problem, a nonlinear factor as shown in Equation (10) was introduced to help describe the transitional change of the deformation modulus.

$$f(t)=\frac{a{t}^b}{1+a{t}^b}$$

(10)

where a and b are fitting parameters which should be positive. Argument t is time. When $t\to 0$, $f\left(t\right)\to 0$. When $t\to \infty$,$f\left(t\right)\to 1$. f(t) is a continuous increasing function which increases from 0 to 1 nonlinearly. The changing trend is controlled by parameters a and b. A nonlinearly changing deformation modulus increment was added to the Merchant model as:

$$E\left(t\right)=\frac{E_2}{f\left(t\right)}=\frac{E_2}{\frac{{\mathit{at}}^b}{{1+\mathit{at}}^b}}$$

(11)

where $E_2$ > 0 is the initial value of the nonlinear component.

By adding certain stress σ₀ to the nonlinear modulus (Equation (11)), the extra creep strain of the nonlinear component can be obtained as:

$$\epsilon ’\left(t\right)=\frac{{\sigma }_0}{E_2}\frac{{\mathit{at}}^b}{{1+\mathit{at}}^b}$$

(12)

By adding Equation (12) to Equation (9), the modified Merchant model can be obtained as:

$$\epsilon \left(t\right)=\frac{{\sigma }_0}{E_0}+\frac{{\sigma }_0}{E_1}\left[{1-e}^{\left(-\frac{E_1}{{\eta }_1}t\right)}\right]+\frac{{\sigma }_0}{E_2}\frac{{\mathit{at}}^b}{{1+\mathit{at}}^b}$$

(13)

Through fitting analysis for all test data, corresponding parameters of the modified Merchant model can be obtained and used for statistical analysis. As shown in Figure 7, fitting results for a group of different loading level (100 kPa) and several repetitive test results were given as an example. Fitting results (means and standard deviations) for all loading levels were shown in

Table 3. Parameter fitting values and statistical characteristics of the improved Merchant model under 100 kPa load.

Creep Parameters	100 kPa		200 kPa		400 kPa		800 kPa		1600 kPa
	Mean	S.D.	Mean	S.D.	Mean	S.D.	Mean	S.D.	Mean	S.D.
E0/MPa	48.738	2.214	41.851	1.614	47.759	0.814	56.812	1.479	67.503	1.325
E1/MPa	910.126	13.808	4027.298	92.566	1120.534	85.481	2647.032	123.936	1070.480	64.753
η1/MPa.h	2332.207	288.536	15,376.760	597.842	13,389.650	1113.850	22,900.400	1828.887	54.018	1.227
E2/MPa	37.705	2.204	60.803	1.014	126.803	11.204	171.750	7.748	60.470	1.450
a	0.517	0.032	0.387	0.008	1.214	0.092	0.934	0.023	0.267	0.008
b	0.200	0.015	0.222	0.017	0.227	0.015	0.230	0.010	0.189	0.009

.

The fitting results demonstrate the superior performance of the modified Merchant model and show the significant variability of the model parameters. Based on the Kolmogorov–Smirnov test, a normal distribution model was chosen to describe the randomness of the parameters. In the modified Merchant model, $E_0$ controls the immediate settlement of soil structure. $E_1$, ${\eta }_1$, ${{\rm E}}_2$, a and b are the parameters that control the damping and stabilizing of the creep deformation. These six parameters are corelated and therefore were treated as a six-dimensional random model as shown in Equation (14):

$$X=\left[E_0 E_1 N_1 E_2 AB\right]$$

(14)

Here, $x_{\left(i\right)}=\left[e_{0i}{ e}_{1i}{ n}_{1i}{ e}_{2i}{ a}_i{ b}_i\right] (i=1,\cdot \cdot \cdot ,n)$ was denoted as a random sample from the statistical population $X\sim N_p(\mu ,\sum ) (p=6)$, and ${X=(e}_{0i}{ e}_{1i}{ n}_{1i}{ e}_{2i}{ a}_i{ b}_i{)}_{n\times p}$ was denoted as all samples from the test results. Joint distribution function of X is

$$F(e_0, e_1, n_1, e_2, a, b)=P\left\{E_0\le e_0, E_1\le e_1, \dots , B\le b\right\}$$

(15)

Joint density function for X is

$$f(x)=\frac{1}{{(\sqrt{2\pi })}^6{\left|\sum \right|}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}\exp \left[-\frac{1}{2}{\left(x-\mu \right)}^{\prime }{\sum }^{-1}(x-\mu )\right]$$

(16)

The mean value vector and covariance matrix of random vector X are

$$\mu =E(X)={\left[E(E_0)E\left(E_1) E\right(N_1)E(E_2)E(A)E(B)\right]}^{{\rm T}}$$

(17)

$$\sum =E\left[\left(X-E\left(X\right)\right){\left(X-E\left(X\right)\right)}^{{\rm T}}\right]$$

(18)

respectively.

Through unbiased estimation method, mean value and covariance of the statistical population were approximated by those of the samples. MATLAB R2020a was used to calculate the mean value vector $\overline{X}$, scatter matrix A and covariance matrix S. The built-in function mvnrnd was used to generate six hundred groups of model parameters for demonstration. Using the acquired parameters, six hundred groups of creep strain at different times can be calculated. After a statistical analysis of the calculated creep strain, the distribution property of the creep strain at different times can be obtained as shown in Figure 8. The calculated creep strain also obeys the normal distribution model. The distribution property for different load levels at a certain time was shown in

Table 4. Distribution property for different load levels at 96 h.

Load/kPa	Time/h	Mean/%	Standard Deviation/%	Strain Range/%
100	96	0.377	0.007	0.36~0.40
200	96	0.653	0.019	0.60~0.72
400	96	1.120	0.035	1.02~1.23
800	96	1.778	0.040	1.63~1.94
1600	96	3.550	0.042	3.43~3.69

.

6. Non-Deterministic Predication Method for Post Construction Settlement of Loess High-Fill Embankment

6.1. Structure of the Proposed Method

Randomness of the creep behavior of compacted loess and other kinds of soils is inevitable. Therefore, it is necessary to consider such randomness in the long-term deformation prediction process for geo-structures such as soil embankments. Based on the modified Merchant model described above, a non-deterministic predication method for post-construction settlement of a loess fill embankment was proposed. As illustrated in Figure 9, the non-deterministic predication method starts with a group of repetitive creep tests that was designed to acquire the random feature of the creep behavior. Then, parameters of the modified Merchant model fitting curve for the test data were obtained through statistical analysis. In this paper, it was assumed that the parameters obey a six-variable normal distribution, and a built-in function mvnrnd of MATLAB was used to randomly generate separate groups of creep parameter values. The time–strain curve for each parameter combination was calculated by substituting the parameter to the modified Merchant model. Then, based on a layer-wise summation method, by proportionally scaling the time–strain curve according to the thickness of each layer and accumulating the deformation of each layer, post-construction settlement of the filling embankment can be calculated. With a group of deformation curves acquired, statistical analysis was then applied to investigate the randomness of the post-construction settlement. The probability contour of the prediction value and the confidence interval envelope under different confidence levels can be obtained.

6.2. Engineering Application

To demonstrate the proposed method in detail as well as prove its validity, an engineering application was conducted to the background engineering case (described in Section 2). Using the modified merchant model and fitting results from Section 5, and followed the proposed non-deterministic predication procedure, probabilistic results of the post-construction settlement were obtained and analyzed.

6.2.1. Calculation Model

The filling height of the background engineering is 61.5 m. Filling construction started in July 2009. Then, after a 112 d downtime period, the construction continued in March 2010. The post-construction monitor started in May 2010. The average bulk density of the filling soil is 20 kN/m³. The fitting result of the modified Merchant model obtained from the laboratory test reveals the relationship between creep strain and time under a certain pressure. For practical usage, a layer-wise summation method was incorporated to realize the layered filling construction. The final post construction creep settlement calculation formula can be expressed as:

$$\begin{array}{c}S\left(t\right)={\displaystyle {\displaystyle \sum }_{i=1}^n}{\mathsf{\Delta }S}_i\hfill \\ \hfill ={\displaystyle {\displaystyle \sum }_{i=1}^n}{\epsilon }_i\left(t\right)H_i={\displaystyle {\displaystyle \sum }_{i=1}^n}\left[\frac{{\sigma }_i}{E_{i0}}+\frac{{\sigma }_i}{E_{i1}}\left({1-e}^{-\frac{E_{i1}}{{\eta }_{i1}}t}\right)+\frac{{\sigma }_i}{E_{i2}}\frac{a_i{t}^{b_i}}{{1+a}_i{t}^{b_i}}\right]H_i\end{array}$$

(19)

where S is the post-construction settlement at a given time t. i is the soil layer number. n is the total soil layer number. $H_i$ is the height of layer i. ε_i(t) is the strain of layer i at time t.

During the calculation, the total 61.5 m height was divided into thirty sublayers with a thickness of 2–3 m according to the actual filling construction. The loading was calculated in accordance with the filling height and bulk density of the filling soil. The settlement at 1d was considered as instantaneous settlement, and the subsequent settlement was considered post-construction settlement.

For a given soil layer, the overlying load was firstly calculated to determine the calculation parameters. Since laboratory tests can only consider finite loading conditions (for example five different loads in this study), it is almost impossible for the load of a soil layer to find matched test settings in practical engineering (for example the load of a soil layer is 150 kPa, but only 100 kPa and 200 kPa were considered in laboratory tests). Therefore, the interpolation method was used to calculate the strain curve of a soil layer for a give load as described by Equation (20):

SC_x = SC_b1 + (x − b1)/(b2 − b1) × (SC_b2 − SC_b1)

(20)

where SC_x is the strain curve of a soil layer under load x. b1 and b2 are the two load boundaries closest to x. SC_b2 and SC_b1 are the calculated strain curve under load b1 and b2, respectively. For example, when the load of a soil layer is 150, its strain curve can be calculated as:

SC₁₀₀ = SC₁₀₀ + (150 − 100)/(200 − 100) × (SC₂₀₀ − SC₁₀₀)

(21)

After calculating the deformation curve for each layer, the total settlement curve was obtained by summing up all layers, and by repeating these steps, multiple settlement curves were acquired, which were used for further statistical analysis.

6.2.2. Prediction Results

Based on the obtained six-dimensional normal distribution model, the creep strain of each soil layer was calculated with Equation (20) with the randomly generated creep model parameters through the built-in function mvnrnd in MATLAB. Figure 10 shows the distribution of obtained six parameters (E₀, E₁, η₁, E₂, a and b) when pressure equals 100 kPa. The deformation curve was then calculated by multiplying the strain curve with layer height. By accumulating all deformation curves of all soil layers, the final deformation curve of the foundation embankment was obtained. By repeating the above steps 600 times with different random parameter sets, the predicted results were obtained and shown in Figure 11, in which 600 independent deformation predictions can be found for each time. In order to verify the validation of the predicted results, three different post-construction prediction methods were applied to the same problem. Two of them used the mean value of the creep test results in this paper as a data source for building a component model [42] and an empirical model [6], respectively. The third method used monitored data for extrapolation [22]. Predicted results of the three methods were shown in Figure 11 as well.

Corresponding monitoring data were also shown in Figure 11. The monitoring point is located 10 m under the surface of the filling embankment. Furthermore, deformation of the original foundation was monitored and excluded so that only the deformation of the filling soil was retained. The red line was placed on the left of the finishing line in Figure 11 because the monitor started several days before the construction finished. It can be seen that the red line which represents the monitored data is located near the center of the predicted result which proved the validity of the proposed method. Additionally, it is clear that the three deformation curves of the comparative methods all fall within the predicted range of the proposed method, which proves the validity of the proposed method on one hand and also illustrates the defect of the existing method that they can only obtain one deterministic prediction despite the fact that the test data are not that certain.

By statistically analyzing the distribution characteristics of the predicted deformations, a probability density function can be fitted at each moment. A contour plot of the probability density of deformation was then obtained and shown in Figure 12. The contour plot graphically illustrates the development of deformation in a non-deterministic way. Along the time axis, the deformation trend was shown, and at a specific moment, distribution characteristics of the deformation can be obtained. To further represent the possible deformation range, different deformation confidence intervals were calculated and shown in Figure 13. The monitored data fall within the confidence interval at a significant level of 0.1.

It can be seen from Figure 11, Figure 12 and Figure 13 that variation of the predicted deformation increases over time and tends to converge in the final stage. When the settlement tends to be stable, the dispersion of settlement also remains unchanged. The engineering prediction range calculated by the non-deterministic prediction method in this paper includes the actual monitoring data, and the calculated prediction value has a large dispersion, indicating that the randomness caused by the geotechnical materials and human factors in construction has a great impact on the engineering deformation, and the randomness should be considered in construction to predict the engineering settlement more accurately.

6.3. Limitations of the Method

It should be noted that the analysis in Section 4 and Section 5 as well as the resulting creep model and the empirical conversion equation is not universal. Different devices, test operators and sample preparation qualities may produce different results even for the same type of soil. These obtained results are only valid for the particular type of compacted loess described in Section 2.1.

The analysis results in Section 6.2 are also limited to background engineering. However, the idea of considering the randomness of the creep behavior and the corresponding non-deterministic prediction method described in Section 6.1 is universal and is the core of this paper. The application scope of the proposed prediction method can be extended to different types of soil or materials that have natural uncertainty by conducting more tests and random analysis on different types of soil.

7. Conclusions

In this paper, uncertainty of the creep behavior of compacted loess was investigated, and a non-deterministic post-construction settlement prediction method was proposed. Firstly, a modified Merchant model was established by introducing a nonlinear element to the original Merchant model to accurately describe the attenuation creep characteristic of compacted loess. Based on the modified Merchant model and a set of repetitive creep tests for compacted loess, a random creep model considering the uncertainty of the creep feature was established, which was defined by a six-dimensional normal distribution model of creep parameters. An empirical transformation equation between staged loading and separated loading mode of the creep test results was also introduced to improve the applicability of the statistical approach. On these bases, a non-deterministic post-construction settlement prediction method was proposed, which was realized by randomly sampling the creep model parameters and statistically analyzing the prediction results calculated by a layer-wise summation method. The proposed method can reflect the uncertainty of the test data in the form of probability contours of the prediction value and the confidence interval envelope under different confidence levels. Finally, the proposed method was applied to a loess filling embankment to predict its post-construction settlement. The measured settlement curve of the engineering falls within the 95% confidence interval of the prediction range, which shows the validity and reliability of the non-deterministic post-construction settlement prediction proposed in this paper. Research in this paper provides a descriptive model for the uncertainty of the creep behavior of compacted loess and a basic approach for the probabilistic design and reliability assessment of filling engineering.