Creep behavior of clayey soil and its model prediction in the Cangzhou land subsidence area

The Cangzhou Plain is an extremely water-deficient area under the dual pressure of population and environment. Because of the scarcity of surface water, the entire region relies mainly on over-exploitation of groundwater to sustain the growing needs of national economic development. Land subsidence poses a threat to the normal operation of municipal facilities such as urban water supply and gas supply, leading to a decrease in the city’s flood control capacity. Land subsidence is mainly due to the compression deformation of soil layers caused by over-exploitation of groundwater, with the compression deformation of clayey soil layers contributing the most to land subsidence. Surface water accounts for about 20% of the total water use in Cangzhou, and the remaining 80% is groundwater¹. With the increase in groundwater extraction, the groundwater level declined and a regional groundwater level depression cone was gradually formed. As a result, the pore water pressure in the aquifer decreased and the clayey soil layers lost water due to compaction, leading to land subsidence. The most severe land subsidence occurred in 2017, with 1400 km² of land experiencing annual subsidence rates exceeding 50 mm².

To mitigate land subsidence, local government had implemented policies since 2005 to limit the extraction of groundwater and to increase the supply of surface water from the Dalangdian Reservoir and the South-to-North Water Diversion Project. As a result, the groundwater level of Cangzhou City has risen significantly, and the land subsidence in the urban area of Cangzhou City has been controlled to a certain extent, but the land subsidence is still in a relatively rapid development stage^2,3. At present, the severe subsidence area in Cangzhou covers an area of about 300 km², mainly distributed in the west and south of Cangzhou. In the central urban area, the maximum cumulative subsidence has surpassed 2.5 m¹.

The relationship between the change of groundwater level and land subsidence in Cangzhou City is shown in Fig. 1. The data were obtained from highly sensitive borehole extensometers provided by the Hebei Cangzhou Groundwater and Land Subsidence National Observation and Research Station².The groundwater level fluctuates and rises after November 2019, but the compression deformation of the aquifer continues to increase, indicating that there is delayed creep in clayey soils. Based on the data of borehole extensometers and groundwater level monitoring, some researchers have suggested the presence of inelastic creep in clayey soils^4,5,6. Although the creep deformation of the aquitard is relatively small compared with the primary consolidation deformation, the land subsidence caused by creep deformation still needs attention due to the large cumulative thickness of aquitards in Cangzhou. Therefore, it is of great significance to analyze the creep characteristics of clayey soil and establish the corresponding creep model for the prevention and control of land subsidence.

Time series of groundwater level depth and cumulative deformation in a confined aquifer in Cangzhou.

Full size image

The macroscopic mechanical behavior of clay is usually represented by constitutive models. Commonly used constitutive models include elastoplastic constitutive models, critical state constitutive models and creep constitutive models^7,8,9. At present, various creep models of soft soil have been developed, including empirical, semi-empirical and semi-theoretical and theoretical models. The theoretical models were generally obtained through rigorous mechanical analysis and derivation, combined with the verification of the results of creep tests. Yuan et al. proposed a generalized elastoplastic-visco model, MIT-SR, which is capable of describing the time-dependent characteristics of creep to shear behavior in clays¹⁰. Ding et al. introduced a constitutive model that can reproduce the shear creep characteristics under complex stress conditions, and the predicted curves of the model show good agreement with the experimental data¹¹. Chen et al. suggested a constitutive model to describe the time-dependent stress–strain behavior of clays at different temperatures, where the viscoplastic strain rate of the soil can be described using three state variables: effective stress, strain, and temperature¹². Such models were rarely applied in practical engineering due to the difficulty in solving them.

Empirical models are widely used in practical engineering applications due to the simplicity of expression, easy access to parameters, etc. Singh and Mitchell proposed an exponential function to describe the stress–strain relationship and a power function to describe the strain–time relationship, on the basis of which a most commonly used three-parameter creep model was developed¹³. Subsequently, the researchers proposed some improved empirical creep models and shear creep theories based on creep test results^{14,15,16,17,18,19,20}. Such models usually use hyperbolic functions and power functions to describe the stress–strain relationship and strain–time relationship. Semi-empirical and semi-theoretical models usually combine the creep empirical relation with mechanical theory to describe the creep behavior of soft soil. Liu et al. established a four-element fractional-order creep model for soft clay based on the fractional-order Burgers creep model, which can more accurately describe the creep curves of soft clay compared to the Burgers model²¹. Chen et al. proposed a four-parameter hyperbolic equation for the relationship between creep volumetric strain, stress, and time, and then a new two-parameter hyperbolic equation for creep was proposed by introducing a reference time, and a five-parameter soft clay creep model was further developed^22,23. The mechanism of such models is clear, but the model parameters are more numerous and their physical significance is unclear or their values are difficult to determine.

Most of the rheological element models belong to the semi-empirical-semi-theoretical models, which have received attention from researchers for describing the complex mechanical behavior of clayey soils using combinations of basic components. Based on the experimental results, many researchers have analyzed stress–strain curves and stress–strain–time curves and proposed improved nonlinear creep constitutive models using combinations of basic elements^{24,25,26,27,28,29,30}. Classical rheological element models include Kelvin model, Maxwell model, Burgers model, etc. As the study continued, the researcher found that the creep properties of some clayey soils could not be described by the above elemental models. The classical element models were improved by connecting the elements in series and parallel to describe the creep properties of clayey soils. For example, the models were improved by introducing nonlinear Hooke element and nonlinear dashpot element^24,25. In recent years, in order to describe the complex creep behaviors of soft clay such as steady-state creep, unsteady-state creep and accelerated creep, the fractional-order nonlinear element had been introduced into the classical rheological element models, whose calculation results were in good agreement with the experimental data^30,31. In addition, numerical and analytical methods were utilized to solve the established constitutive models^32,33, which provided a reference for the calculation and analysis of deformation of multi-layer clay soil. In general, the more creep properties are reflected by a creep model, the more parameters the model needs, which often cannot be directly determined by conventional laboratory or field tests. In order to reduce the testing cost and improve the accuracy of model parameter identification, some evolutionary algorithms have been applied in model parameter identification and optimization and have achieved good results^{34,35,36,37,38,39}.

The research objective of this paper is to develop a nonlinear model reflecting the creep properties of deep clayey soils. Firstly, creep tests were carried out on clayey soils in aquitards in a typical subsidence zone in Cangzhou. Secondly, on the basis of analyzing the creep and deformation properties of clayey soils, a nonlinear creep model of NCE_CS is proposed by improving the Kelvin body in the Merchant model and replacing the elastic element by the power function element, and optimizing and solving the model parameters using genetic algorithm. This work will provide a theoretical basis for the quantitative calculation of creep deformation of aquitards in Cangzhou area.

Study area

Cangzhou City is located in the eastern part of the North China Plain (Fig. 2a) and has a flat topography. The predominant stratigraphy in this area is loose unconsolidated Quaternary sediments with the thickness from 200 m to over 600 m (Fig. 2b, c), consisting primarily of sand, silt, and clay. The predominant type of groundwater is pore water in the loose Quaternary sediments, and the aquifer systems typically have multiple sandy aquifers separated by extensive and thick silt or clay aquitards (Fig. 2c). From top to bottom, there are mainly four aquifer groups, which roughly correspond to the Holocene (Q₄), the Upper Pleistocene (Q₃), the Middle Pleistocene (Q₂) and the Lower Pleistocene (Q₁), respectively. The bottom depths of the first, second, third and fourth aquifer groups are 40~60 m, 120~170 m, 250~350 m, and 350~550 m, respectively. Groundwater in the first and the second aquifer groups has unconfined or semi-confined hydraulic characteristics, and groundwater in the third and the fourth aquifer groups is confined.

Study area and the location of sampling boreholes, and a typical hydrogeological cross section A–A′.

Full size image

Geotechnical testing

The thickness of the Quaternary clay layer in Cangzhou City is large, reaching a thickness of nearly 200 m. Since the shallow groundwater in Cangzhou is dominated by saline/brackish water, the aquifer groups below the first aquifer group are the main groundwater abstraction layers, and their compression is also the reason for the development of land subsidence. In this paper, four undisturbed clayey soil samples below 65.7 m were selected for the geotechnical creep tests, which were taken from the urban area of Cangzhou, and the sampling location is shown in Fig. 2c.

The characteristics of the four undisturbed samples (Fig. 3) are as follows: (1) The clay is mainly yellowish-brown and yellow–brown in color, and its composition is nonuniform, containing small sand grains and white calcium sediments; (2) The liquid index is less than 0.25, and it is in the state of hard-plasticity or hardness. The basic physical parameters are shown in Table 1.

Samples of undisturbed soil used in the tests. (a) N1 (b) N2 (c) N3 (d) N4.

Full size image

Under groundwater abstraction conditions, the groundwater level continues to fall, forming multiple localized cones of groundwater depression, which gradually connect to form a regional depression cone. The groundwater level depth in the center of the depression cone in Cangzhou area has exceeded 90 m. The changes in the stress states of the deep clayey soil due to the drop of groundwater level are shown in Fig. 4.

Schematic diagram of the stress state change in the soil due to groundwater level decline.

Full size image

For the deep clayey soil units, the additional stress caused by drop in groundwater level can be expressed as⁴⁰

where Δσ_z is the vertical additional stress; γ_sat, γ′ and γ_w are the saturated unit weight and buoyant unit weight of soil and unit weight of water, respectively; Δh is the value of groundwater level decline.

Clayey soil units will have a tendency to deform laterally when subjected to vertical additional stress Δσ_z, at which point the horizontal pressure changes from static earth pressure σ₀ to passive earth pressure σ_p. According to Rankine’s earth pressure theory (Eq. 2), the passive earth pressure can be expressed as⁴⁰

where, γ , γ_sat and γ′ are the natural unit weight, saturated unit weight and buoyant unit weight of soil, respectively; k_p is the passive earth pressure coefficient, greater than 1; z is the depth of the soil element; c and ϕ are the cohesive force and internal friction angle, respectively.

When the depth z is large, the additional stresses Δσ_z induced by the drop of groundwater level are much smaller than σ_p, so the lateral deformation can be neglected. In this way, it is reasonable to use one-dimensional compression test to study the creep behavior of deep clayey soil.

The one-dimensional creep tests were carried out with reference to the Standard for Geotechnical Test Methods (China Planning Press, 2019), and the test apparatus was a WG-type high-pressure consolidation apparatus produced by Nanjing Soil Instrument Factory. The specimen was 2.0 cm high and had a cross-sectional area of 30.0 cm², and the specimen was drained on both sides during the consolidation process. The specimens were loaded in a stepwise loading mode, and the loading scheme of specimens N1, N2 and N4 was as follows: 100 kPa → 200 kPa → 400 kPa → 800 kPa → 1600 kPa → 3200 kPa, and the loading time for each level was 4d, which lasted for a total of 24d. The loading scheme of specimen N3 was as follows: 100 kPa → 400 kPa → 800 kPa → 1600 kPa → 3200 kPa, and the loading time for each level was 8d, which lasted for 40d in total. The laboratory temperature was controlled at 24 ± 1 °C during the test to ensure that the outside temperature had minimal effect on the test results.

Analysis of creep curves

Figure 5 shows the strain–time curves of N1~N4 specimens under each level of load. As shown in Fig. 5a, the strain–time curves under the whole loading present a stepped shape, which has an obvious turning point at the initial stage with the increase of each loading, and then enters the steady-state creep stage. The grading loading testing (Fig. 5b) shows that the total strain (varepsilon) of each loading is composed of three parts: instantaneous strain (varepsilon_{i}), primary consolidation strain ({varepsilon }_{p}) and creep strain (varepsilon_{c}) (Eq. 3):

Creep testing curves of the clayey soil samples N1–N4.

Full size image

The reasons for the turning point of the strain–time curve are as follows. The specimen generates excess pore water pressure under loading, and the effective stress increases during the gradual dissipation of pore pressure, leading to the primary consolidation deformation which accounts for a large proportion of the total deformation. After the pore water is completely discharged, the deformation is mainly dominated by the creep of soil particles and cohesive water, accounting for a small proportion of the total deformation. According to Casagrande’s method, e-log(t) curves were made to divide the primary and secondary consolidation. Casagrande believed that the end time of primary consolidation in clayey soils is about 200 min. The test results for specimens N1~N4 showed that the end time of primary consolidation was about 400 min, and thus 400 min was chosen as the time of creep occurrence.

For specimens N1~N4, because the duration of stepwise loading is the same and the deformation is in a steady development trend, so it is more appropriate to use the Boltzmann superposition principle to obtain the creep curve under different loads. The creep curves under each loading was obtained by the coordinate translation method, as shown in Fig. 6. When the loading time t = 400 min, the creep strain ε_c = 0. The creep curves of the clayey soil specimens have obvious nonlinear characteristics, creep deformation increases gradually with the increase of time, but the increase rate generally shows a decreasing trend. It can also be seen from Fig. 6 that the amount of creep increases gradually as the load increases.

ε_c–t relations in the creep stage under grading loading. (a) N1 (b) N2 (c) N3 (d) N4.

Full size image

Analysis of isochronous stress–strain curves

The isochronous stress–strain curves of soil samples N1~N4 were obtained by using the Boltzmann superposition principle, as shown in Fig. 7, the curves under different loading have similar shapes, and the strain increases with the increase of time for the same load, but the isochronous curves gradually get close to form a curve cluster. The creep deformation increases with the increase of stress, and the stress–strain shows a nonlinear relationship, which should be taken into account when constructing the creep models.

Isochronous stress–strain curves of the soil samples N1~N4. (a) N1 (b) N2 (c) N3 (d) N4.

Full size image

Classical rheological models

The basic elements commonly used to describe the rheology of materials include Hooke elastomer, Newton viscous body and St. Venant plastic body. By combining these basic elements, the various creep models can be obtained, which can be used to present the relationship between stress, strain and time of rock and soil mass. The classical models and mathematical expressions are shown in Table 2. Because the basic components of these models have clear physical meanings, which is convenient for the in-depth analysis of the research object, this type of models are more widely used in the creep study of rock and soil mass.

An improved nonlinear creep model

The classical creep models above have the following shortcomings: (1) They are not effectively combined with the consolidation theory of soil mechanics, but based on the viscosity, elasticity and plasticity of general materials, which cannot reflect the deformation caused by consolidation and drainage of clayey soil, and the results calculated by the above models are inconsistent with the actual observation results; (2) They cannot reflect the three stages of instantaneous deformation, primary consolidation and secondary consolidation (creep) in the consolidation theory; (3) When t = 0, the creep strain is theoretically 0, however, the results calculated by the above models all have transient strains, such as the Burgers model, which calculates the creep strain at t = 0 to be (varepsilon = {sigma mathord{left/ {vphantom {sigma {E_{H} }}} right. kern-0pt} {E_{H} }}), which is inconsistent with reality. Considering the creep behavior of clayey soil in land subsidence areas, the NCE_CS nonlinear creep model is proposed by replacing elastic elements in Merchant model with power function elements and improving Kelvin body in Merchant model to consider the beginning time of creep in consolidation theory.

Improvement of Kelvin model

Kelvin model consists of a Hooke elastomer and a Newton viscous body in parallel (Table 2). When time t → ∞, the strain tends to be a constant, reflecting the property of steady-state creep of the soil. The strain after unloading takes some time to decrease to 0, reflecting the nature of elastic aftereffects. The Kelvin model can be expressed as²⁶

where ({varepsilon }_{kc}) is the creep of clayey soil under load (sigma_{K}), E_H is Hooke elastic modulus, η is the viscosity coefficient in MPa·min, and (dot{varepsilon }_{kc}) is the creep that occurs per unit time and reflects the speed of creep deformation.

In the strain–time natural logarithm coordinate system, the creep rate can be expressed as:

where t is any time in the creep phase; t₀ is the start time of creep stage, here 400 min is taken; min is the time unit minute.

Substituting Eq. 5 into Eq. 4:

Equation 7 is derived by transposition of Eq. 6:

When t = t₀, creep begins, and let (varepsilon_{kc} left( {t_{0} } right) = 0); By integrating both sides of Eq. 6 in the range of t₀ to t, Eq. 8 can be obtained:

The simplification of Eq. 8 is as follows:

By sorting out Eq. 9, the improved Kelvin equation is obtained:

The improved Kelvin creep model is more concise and is able to take into account the time of occurrence of the secondary consolidation (creep) phase in consolidation theory, making the division between creep and primary consolidation deformation clearer.

Power function element

The creep of clayey soil can be highly fitted by a power function¹⁵. The proposed power function element, as shown in Fig. 8, is similar to a viscoelastic element and can reflect the nonlinear behavior of creep. Its equation can be expressed as Eq. 11.

Power function element.

Full size image

Equation 11 is expressed in stress form as Eq. 12:

where (varepsilon_{pc}) is the calculated creep, (sigma) is the magnitude of the load, and a, b, c are model parameters.

Using this power function element, the creep at the moment t = t₀ can be calculated to be 0, and the effects of variable dimension are eliminated by using t/t₀ as the creep time.

NCE_CS model

Combined with the one-dimensional creep test results of clayey soil in the land subsidence area of Cangzhou, the power function element replaces the Hooker’s elastic element in the Merchant model and is connected in series with the improved Kelvin body to obtain the nonlinear creep NCE_CS model, which is schematically shown in the model structure in Fig. 9.

NCE_CS nonlinear creep model.

Full size image

The NCE_CS model is expressed as Eq. 13:

Considering that the isochronous stress–strain curves (Fig. 7) have nonlinear characteristics, the modulus (E_{H}) in Eqs. 10 and 12 is different from the elastic modulus (E_{se}) in Eq. 13, which varies with the change of loading and can be calculated from the segmented compression modulus on the e-log(P) curves.

where (varepsilon_{c} left( t right)) is the creep under the stress (sigma); a, b, c are model parameters;(E_{se}) is the segmental modulus at each level of loading; (E_{si}) is the segmental compression modulus under load level i; (sigma_{i + 1}) and (sigma_{i}) are the magnitude of the load of the (i + 1)th and i-th level, respectively; (e_{i + 1}) and (e_{i}) are the pore ratio under the load of the (i + 1)th and i-th level, respectively; (e_{0}) is the initial pore ratio; μ is the Poisson’s ratio.

The mechanical behavior of clayey soils exhibits strain rate sensitivity, which means that their strength and deformation vary with different loading rates. Elastoplastic constitutive models usually introduce strain rate effects indirectly through over-stress theory, thus limiting their applicability. Unlike elastoplastic models^41,42, the proposed model not only captures the time-dependent deformation under long-term loads, but also directly reflects the strain rate sensitivity through the viscous Kelvin element.

Determination of model parameters

The parameter solution of NCE_CS model is actually a multivariate function parameterization problem, and the genetic algorithm can find the global solution more accurately, which is very suitable for the optimization of multivariate function parameterization. The NCE_CS model includes four parameters (a, b, c, η), which can be solved by the genetic algorithm through the module sko.GA, which is written in Python language. The algorithm is a simulation of the biological evolution process, taking creep (varepsilon_{c} left( t right)) as the objective function (Eq. 13), time t and stress σ as variables, and optimizing the parameters by random search.

The solution steps are as follows (Fig. 10): (1) Initialize the population; the population number is set to 5000 in this paper. (2) The individuals are assessed to determine their fitness. (3) Selection and crossover; the individuals are selected according to the fitness, then the selected individuals are crossed and new individuals are produced. (4) Mutation; the diversity of samples is increased through variation, and the mutation probability is set at 0.001 in this paper. (5) Repeat the above steps; the calculation is stopped when the maximum number of iterations (1000 times in this paper) is reached or the optimal solution is obtained. The modulus E_se is obtained by the consolidation test data and Eqs. 14 and 15. The calculation results of the model parameters are shown in Table 3.

Flow chart of genetic algorithm solution.

Full size image

The parameters in Table 3 are substituted into Eq. 13, and the creep calculation results and SR under various loads are obtained. (SR = frac{{varepsilon_{c – cal} left( t right)}}{{varepsilon_{c – tes} left( t right)}}) is a ratio between the calculated value (varepsilon_{c – cal} left( t right)) and the testing value (varepsilon_{c – tes} left( t right)) at any time point, and the ratio can reflect the accuracy of the calculated results. The variation of SR values of specimens N1~N4 with time under different loads is shown in Fig. 11, from which it can be seen that the SR values are generally large at the initial moment and when the load is small, and the remaining SR values varying around 1.0. The maximum value of SR is 1.27 and the minimum value is 0.77, 85.1% of the SR values are between 0.8 and 1.2, which indicates that the model predicts well.

Changes of SR values with time under different loads. (a) N1 (b) N2 (c) N3 (d) N4.

Full size image

Dependence of parameters a, b and c on consolidation pressure

In order to make the model more widely applicable, this paper tries to find the relationship between the model parameters and other physico-mechanical indexes. Imitating the definition of the overconsolidation ratio, the load action ratio W is defined as the ratio of the soil’s self-weight stress to the total pressure currently borne by the soil. The relationship between parameters a and b and the load action ratio W is shown in Fig. 12, approximately in the form of power function, indicating that the parameters a and b decrease with the increase of W. When W is between 0 and 4.0, the rate of decrease of parameters a and b is large, and the rate after W > 4 is small. The nonlinear curve fitting function in Origin software was utilized to fit the experimental data, and the relationship between parameters a, b and W is shown in Eqs. 16 and 17, whose coefficients of determination are 0.855 and 0.955, respectively, indicating a good fitting effect.

Relations between parameters a, b and W.

Full size image

To make the variables dimensionless, the load (sigma) divided by standard atmospheric pressure ({text{P}}_{{text{a}}}) ((LR = {raise0.7exhbox{$sigma $} !mathord{left/ {vphantom {sigma {{text{P}}_{{text{a}}} }}}right.kern-0pt} !lower0.7exhbox{${{text{P}}_{{text{a}}} }$}})) is defined as the load-to-pressure ratio. The parameter c decreases with the increase of LR, which is approximately in the form of a power function (Fig. 13). The relationship between parameter c and LR can be expressed in Eq. 18 with a coefficient of determination of 0.968.

Relations between parameters c and LR.

Full size image

The influence of parameter (eta) on the calculation results

The analysis shows that the regularity of the relationship between model parameters such as viscosity and load-action ratio and physical property indices is poor. It is necessary to analyze the effect of (eta) on calculation results of the creep. The values of (eta) in Table 3 range from 2.120 to 2.979, with a mean value of 2.6. The value of (eta) can be taken as 2.0, 2.2, 2.4, 2.8, 3.0 and 3.2, respectively, and the creep is calculated according to Eq. 13, and then compared with the creep calculated at (eta) = 2.6. The error is calculated according to Eq. 19.

where (P_{er}) is the average of the creep calculation errors for different values of (eta). N is total number of comparisons;(varepsilon_{eta ,i}) is the creep calculated by different (eta) at the i-th comparison; (varepsilon_{2.6,i}) is the creep calculated with (eta = 2.6) at the i-th comparison.

As the deviation from (eta = 2.6) increases, the error (P_{er}) of creep calculation also increases, up to a maximum of only 2.46% (Fig. 14). When (eta = 2.4) and (eta = 2.8), the error is close to 0. This indicates that the variation of (eta) in the range of 2.0~3.2 has little effect on the creep calculation results, and so let (eta = 2.6).

Errors in creep calculation (Eq. 19) when (eta) takes different values.

Full size image

Equations 16, 17, 18 and (eta = 2.6) are substituted into Eq. 13 to obtain the expression of the nonlinear creep NCE_CS model (Eq. 20):

Verification and application of the model

Comparison of predicted results with testing data

In order to verify the reliability of the NCE_CS model, the calculated results of the NCE_CS model were compared with the test data and the calculated results of the Kelvin model, Maxwell model, Merchant model, and Burgers model, respectively, in terms of the creep curve and the isochronous stress–strain curve. The comparison results of specimen N3 are shown in Figs. 15 and 16. Based on the experimental data, genetic algorithm was used to find the parameters of the above four classical creep models, as shown in Table 4.

Creep curves calculated by different creep models versus test data.

Full size image

Stress–strain curves calculated by NCE_CS model, Merchant model versus test data.

Full size image

Figure 15 shows that the creep curve of the NCE_CS model is more consistent with the test data, and the creep shows a nonlinear steady state development trend. The calculation results of the Merchant model are obviously smaller than the test data, and the creep growth is faster at the beginning, and then increases slowly at the later stage of the creep (Fig. 15c). The calculation results of the Kelvin model are very different from the test data, and the results evolve into a nearly horizontal straight line with the increase of time (Fig. 15a). The calculation results of Maxwell model and Burgers model have a big difference with the test data, and the calculation results show a linear increase with the increase of time, and the slope increases with the increase of load (Fig. 15b, d). Since the creep curves calculated by the Kelvin model, Maxwell model, and Burgers model differ significantly from the experimental data, only the isochronous stress–strain curves of the NCE_CS model and the Merchant model are compared with the experimental data. As shown in Fig. 16, there is a discrepancy between the results calculated by the Merchant model and the experimental data, and the shapes of the curves are also inconsistent. The creep isochronous stress–strain curves calculated by NCE_CS model are in good agreement with the experimental data, which indicates that the NCE_CS model is able to describe the creep characteristics of clayey soil in Cangzhou.

Application of NCE_CS model

In order to investigate the causes of land subsidence, hydrogeological drilling and sampling were conducted in Renqiu City in northwest of Cangzhou (Fig. 2b). The drilling reached a depth of 400.0 m (Fig. 17a), in which clayey soil samples were collected. The basic geotechnical test and compression creep test were carried out on the clayey soil sample R1 buried at a depth of 58.86 m. The physical indexes of the soil sample R1 are shown in Table 5. The compression creep tests were carried out in a MIS-232 high pressure consolidation instrument from Japan with applied loads of 0.325 MPa, 0.65 MPa, 0.975 MPa, and 1.3 MPa, and the loading duration was 112 h. The stepwise loading strain–time curves (Fig. 17b) of specimen R1 are similar to those of specimens N1~N4, and the creep deformation under different loads shows a steady state development trend as time increases.

Borehole histogram and stepwise loading creep curves in Renqiu.

Full size image

By calculating and analyzing the test data, the values of E_se, W and LR at all levels of loading were obtained (Table 6). The value of (eta) is still taken as 2.6.

Comparing the clayey soils of Renqiu City and Cangzhou City, due to differences in composition, deposition environment, and groundwater level, the power function relationship between parameters a, b, c and W, LR in Eq. 20 remains unchanged, but the values of the constant terms in Eq. 16, 17, and 18 will be changed. The parameters of the NCE_CS model were solved using genetic algorithm, and Eq. 21 was finally obtained.

The creep curves (Fig. 18a) at all levels of loading were calculated using Eq. 21, which are overall consistent with the test data, with a large deviation only in the initial stage at a load of 0.325 MPa. Figure 18b shows that the maximum value of SR, which reflects the calculation accuracy, is at the initial stage of creep under the load of 0.325 MPa. The maximum value of SR is 7.694, the minimum value is 0.957, and the percentage of SR values between 0.9 and 1.2 is 78.95%, which indicates that the NCE_CS model can reflect the creep law of the clayey soil in this area.

Calculation results of NCE_CS model versus test data of clayey soil in Renqiu.

Full size image

This paper analyzes the characteristics of groundwater level change and the physical state and stress state of deep clayey soil in the land subsidence area of Cangzhou, and investigates the creep behavior of deep clayey soil by using compression tests. Combined with the test data and theoretical analysis, the study proposes a creep model of clayey soil and draws the following conclusions.

1. 1.

   The test results show that the total strain under stepwise loading consists of three parts: instantaneous strain, primary consolidation strain and creep strain. According to the Boltzmann superposition principle, the creep curves and isochronous stress–strain curves under stepwise loading are obtained, and both of them have obvious nonlinear characteristics, and the creep curves as a whole show a steady state development trend.

2. 2.

   Based on the Merchant model, the NCE_CS nonlinear creep model was established by replacing the elastic element in the model with a power function element and improving the Kelvin body of the Merchant model so that the onset time of creep can be considered. The NCE_CS model introduces the onset time of the creep stage, which makes the division of the primary consolidation and creep stages clear and effectively reflects the nonlinear creep characteristics of clayey soils.

3. 3.

   The NCE_CS model has five parameters, among which the modulus E_se can be obtained by one-dimensional compression test, and the remaining four parameters a, b, c, and η can be calculated by genetic algorithm. The calculated results of the NCE_CS model were compared with those of the test data and four classical creep models. The results show that, compared with the classical creep models, the creep curves and isochronous stress–strain curves calculated by the NCE_CS model are in better agreement with the test data, indicating that the NCE_CS model can better reflect the creep deformation characteristics of clayey soils.

4. 4.

   To further verify the applicability of the NCE_CS model, the creep curves of the clayey soils in another subsidence area, located in Renqiu City in northwest of Cangzhou, were obtained by using this model, with most of the SR values ranging from 0.9 to 1.2, indicating that the model calculations are highly accurate. Once again, it shows that the NCE_CS model can effectively reflect the creep properties of clayey soils in Cangzhou area.