Finite Difference Method for the One-Dimensional Non-linear Consolidation of Soft Ground Under Uniform Load

Initial stress and additional effective stress distributions in soil greatly influence the degree of ground consolidation when calculating one-dimensional soft clay ground consolidation in deep soil. The one-dimensional non-linear consolidation governing equations of soft ground under uniform load are derived and solved with the finite difference method. This method is based on the assumptions that the initial stress in soil varies with the ground depth and that the additional effective stress caused by external loads changes with both the ground depth and consolidation time and the hyperbolic model of the soil stress–strain relationship. Formulas for the degree of consolidation and the settlement of the ground are presented. A case study shows that the degree of consolidation in the ground calculated with the finite difference method agrees well with the traditional analytical solution, and the computational efficiency of the finite difference method can be effectively improved when the segmental calculation method is used throughout the consolidation process. The results of another example show that the settlement of the ground calculated with the finite difference method agrees with the in situ data. The suggested method can greatly simplify the consolidation calculation and has a high application value in engineering.

Initial stress and additional effective stress distributions in soil greatly influence the degree of ground consolidation when calculating one-dimensional soft clay ground consolidation in deep soil. The one-dimensional non-linear consolidation governing equations of soft ground under uniform load are derived and solved with the finite difference method. This method is based on the assumptions that the initial stress in soil varies with the ground depth and that the additional effective stress caused by external loads changes with both the ground depth and consolidation time and the hyperbolic model of the soil stress-strain relationship. Formulas for the degree of consolidation and the settlement of the ground are presented. A case study shows that the degree of consolidation in the ground calculated with the finite difference method agrees well with the traditional analytical solution, and the computational efficiency of the finite difference method can be effectively improved when the segmental calculation method is used throughout the consolidation process. The results of another example show that the settlement of the ground calculated with the finite difference method agrees with the in situ data. The suggested method can greatly simplify the consolidation calculation and has a high application value in engineering.

INTRODUCTION
The traditional Terzaghi one-dimensional consolidation theory assumes that external loads are applied instantaneously and does not take into account that external loads are often linear, graded, or cyclic in practical engineering. Additionally, the non-linear characteristics of the stressstrain relationship of the soil itself are not considered. Therefore, many scholars have carried out improvement studies. Schiffman (1958) first analyzed the problem of the one-dimensional consolidation of a foundation under variable loads. Later, Wilson and Elgohary (1974), Favaretti and Mazzucato (1994), Xie and Pan (1995), and Zhao and Shi (1996) conducted further research on this problem and obtained many useful conclusions. Lai and Song (2005) conducted soil consolidation and drainage tests under special stress paths, and the stress-strain relationship of soil obtained was consistent with the fitting results of the log-curve non-linear model, sinusoidal non-linear model, and traditional hyperbolic model. Yu et al. (2008) used a modified hyperbolic model to fit the quantitative relationship between the secondary consolidation coefficient and the pressure of soil under normal consolidation and proposed a modified method for calculating the secondary consolidation settlement of the normal consolidated soil. He and Jiang (2009) proposed a hyperbolic prediction model that reflects the semicubic nonlinear relationship between the settlement rate and the remaining settlement. In addition, in the study of one-dimensional nonlinear consolidation theory, Davis and Raymond (1965), Barden and Berry (1965), and Xie et al. (2006) considered the physical non-linear characteristics of soil based on the relation elgσ ′ . However, the value of the initial consolidation pressure is crucial to the calculation of the foundation settlement by the e -lgσ ′ model. In view of this, Wei (1987) believed that, if the initial pressure and repressure could not be reasonably distinguished, then the hyperbolic model could be used to calculate the foundation settlement, and the calculation result was better than that of the e -lgσ ′ model. Xu (1987) also experimentally proved that the hyperbolic model could better simulate the constitutive relationship of soft clay. On this basis, Shi et al. (2001) obtained the one-dimensional consolidation theory of the hyperbolic model at the time of instantaneous loading, and the calculated results were in good agreement with the laboratory test. Xie et al. (2010Xie et al. ( , 2016 obtained the analytical solution for the one-dimensional consolidation of clayey soils and double-layered structured soils with a threshold gradient. Lo et al. (2016) presented a one-dimensional consolidation theory in unsaturated soils under cyclic loading. Wang et al. (2017), Liu et al. (2018), Zou et al. (2018), and Zhao et al. (2020), respectively, obtained analytical solutions for the one-dimensional consolidation of clayey soils considering different load and boundary conditions. When considering time-dependent loading, Deng et al. (2019) presented closedform solutions for one-dimensional consolidation in saturated soils, Moradi et al. (2019) analyzed the one-dimensional consolidation of a multilayered unsaturated soil under partially permeable boundary conditions, and Zhou et al. (2020) obtained an analytical solution for a classical one-dimensional thaw consolidation model.
However, due to the small sample used in the test, the applicability of this theory in the calculation of deep soft clay could not be verified. Zhang and Sun (2007) established the one-dimensional consolidation theory of foundations under variable loads based on the hyperbolic model, and Zhou et al. (2013) presented a finite difference model for a onedimensional electro-osmotic consolidation. However, the theory assumed that the initial effective stress and the additional stress of soil remained unchanged along the depth, which was different from the actual situation. Therefore, this paper combines previous research results with the hyperbolic model of the stress and strain of soil as a starting point. This method assumes that the initial stress changes along the depth of the foundation, and at the same time, the load caused by the additional effective stress changes with time and depth. Additionally, the uniformly distributed load is derived via the control equation of the one-dimensional consolidation of soil, and the finite difference method is introduced into the degree of consolidation and settlement calculation formula.

Basic Assumptions
The basic assumptions are as follows: the soft soil is saturated, the thickness of the compressed soil layer is H, the surface of the soft soil is permeable, and the drainage conditions of the bottom surface can be divided into permeable and impermeable. Considering the simulation of vehicle load, the impulse load and simple harmonic load, which varies with time, can be selected. When the vehicle load is transferred to the embankment, the impulse load and simple harmonic load could be uniformed. Besides, the embankment is filled step by step, and the embankment load, which can be simulated to multiple graded loads, also varies with time. All these loads above can be simplified into the uniform load q top (t), which varies with time.
The distribution of the water pressure in the foundation is shown in Figure 1, where γ w (h 1 + Hz) is the hydrostatic pressure at depth z, p w is the total water pressure, and u is the pore water pressure generated by the evenly distributed embankment load q top (t) in the soil. Frontiers in Earth Science | www.frontiersin.org a. The stress-strain relationship of the soil satisfies the hyperbolic model: where σ ′ is the effective stress, E 0 is the initial compression modulus of the soil, ε is the strain, and m is the slope of the σ ′ /ε ∼ σ ′ curve.
b. Let the permeability coefficient k be proportional to the compression coefficient a v : where γ w is the unit weight of water, γ w = ρ w g. The parameter ρ w is the density of water, c v is the consolidation coefficient, e 0 is the initial void ratio, and const is a constant.
c. The initial effective stress of the soil σ 0 ′ (z) changes with depth z. d. Considering the change in the external load with time (as shown in Figure 2A), the additional stress of the soil caused by the external load also changes with time, and the additional stress will change along the depth, so the additional stress can be expressed by q(t,z) (shown in Figure 2B).

Derivation of the Consolidation Governing Equation
Based on basic assumptions a and b, the one-dimensional consolidation governing equation can be obtained as: According to the effective stress principle, the stress expression at depth z is The two sides of Equation (4) are opposite and partial derivatives, and we obtain: Substituting Equations (5) and (6) into Equation (3), the one-dimensional non-linear consolidation governing equation of the foundation under uniformly distributed load is Boundaries and Initial Conditions a. At the initial moment, the upper load is conveyed by the pore water pressure, and its mathematical expression is b. After loading, when the ground surface is a permeable layer, the pore water pressure on the bottom surface is zero.
c. When the ground surface is an impermeable layer after loading, the surface pore water pressure is zero, and its mathematical expression is d. When the top surface of the foundation is permeable and the bottom surface is impermeable, the pore water pressure of the soil is equal to that of the underlying layer.

Difference Method Format
According to the basic principle of the finite difference method, with depth z as the vertical axis and time t as the horizontal axis, a calculation area for foundation consolidation analysis is established. The depth z is divided into M equal divisions, the time t is divided into N equal divisions, and a difference method calculation grid is established. As shown in Figure 3, the grid node represents the pore water pressure u, where u a,b represents the value of the pore water pressure at depth b in the ground at the moment of time a.
The initial and boundary conditions in the difference equation can be expressed as Frontiers in Earth Science | www.frontiersin.org For the governing Equation (7), the Crank-Nicolson format is used to discretize the following difference format: where t = T/N , h = H/M , T is the total consolidation time, and M and N are the number of grids of the depth and consolidation time, respectively. To simplify the calculation, Equation (12) is linearized; then, Linearization will inevitably bring errors. To minimize the errors, the difference grid can be refined. Assume Then, the formula for calculating the difference method under two boundary conditions can be written as follows: When the bottom surface is permeable, A is an (M -1) × (M -1) order matrix, u is an (M -1) × 1 order matrix, and B is an (M -1) × 1 order matrix. When t = a t, then When the bottom surface is impermeable, A is an M × M order matrix, u is an M × 1 order matrix, and B is an M × 1 order matrix. When t = a t, then

Calculation Formula of the Consolidation Degree
According to the uniaxial compression theory, the calculation formula of the ground settlement is The formula for calculating the average consolidation degree of a foundation defined by settlement can be written as In addition, the average consolidation degree of the foundation, as defined by the average pore pressure, can also be expressed as

Theoretical Verification
Examples are cited in the literature (Zhu and Yin, 1998) for verification. The calculation parameters are as follows: the thickness of the compressed soil layer is 5 m, and the bottom surface is a permeable layer; the load time linearly increases from q a = 0 kPa to q a = 300 kPa after 50 days and then remains unchanged; the consolidation coefficient of the soil is c v = 5.7888 × 10 −3 m 2 /day, and the initial compression modulus is E 0 = 1,687 kPa. The analytical solution in Zhu and Yin (1998) assumes that the additional stress caused by the external load is linearly distributed with time and depth, and the stress-strain relationship of the soil is a linear elastic model. Corresponding to the special case of the difference method solution in this article, which is m = 0, Equation (1) can be written as In the following, the finite difference method (hereinafter referred to as the method in this paper) proposed earlier is used to calculate the average consolidation degree of the foundation.
To improve the calculation efficiency, the difference calculation is divided into two sections in the time dimension: the first section is 100 days in length, which is divided into 500 parts, t = 0.2 day; the length of the second period is the remaining 1,800 days, which is divided into 1,800 parts, t = 1 day. The average consolidation degree of the foundation calculated by this method and the analytical solution calculation results in Zhu and Yin (1998) are summarized in Table 1.
From Table 1, it can be seen that the maximum relative error of the calculation results of the method in this paper and the analytical solution in Zhu and Yin (1998) is only 0.0106%, indicating that the method in this paper is reasonable and feasible and can meet the calculation accuracy requirements.

ILLUSTRATIVE EXAMPLES
In order to verify the feasibility of the method for the calculation of the degree of consolidation of the foundation under the conditions of a variable load and one-dimensional non-linear consolidation, a calculation example in Favaretti and Mazzucato (1994) was selected for analysis. This example is a weak foundation project located under the silo of Ca'Mello. The change law of the load applied to the soft foundation with time is shown in Figure 4. Favaretti and Mazzucato (1994) carried out on-site observations of the settlement of the soft foundation and measured the consolidation coefficient of the soil (c v = 0.1296 m 2 /day) and the change law of the soil stress and strain during loading by experiments. In this paper, the compression characteristics of the hyperbolic model of the soil (E 0 = 270 kPa, m = 0.9) are obtained by fitting the soil properties from on-site observations in Favaretti and Mazzucato (1994). The method in this paper is used to calculate the ground settlement, which are compared with the one calculated by an analytical solution in Rahal and Vuez (1998) and the one measured in Favaretti and Mazzucato (1994). The results are shown in Figure 5.
It can be seen from Figure 5 that the fluctuation trend of the settlement curve calculated by the method in this paper is basically consistent with the measured results, while the analytical solution settlement curve in Rahal and Vuez (1998) lags behind the measured results, indicating that the method in this paper is more reasonable and feasible.
In fact, the analytical solutions of the foundation settlement are often derived under the standard load form, which cannot flexibly deal with the problem of variable loads in actual engineering. When using the theories in Rahal and Vuez (1998) to calculate the settlement of the soft foundation, the external load needs to be simplified to a standard simple harmonic load. This simplification will inevitably lead to human error, which will cause the theoretical settlement curve to be out of sync with the measured curve. The method in this paper only needs to input the function expression of the external load into the program. Impulse loads, simple harmonic loads, or multiple graded loads can be accurately simulated, which simplifies the load form conversion step, reduces human error, and greatly simplifies the calculation process of the degree of consolidation and the settlement of the foundation. Besides, the accuracy can be improved by calculating in sections in the time dimension, which has the irreplaceable advantage than the analytical solution.

CONCLUSIONS
a. The presented method can adapt to various variable load situations, can greatly simplify the consolidation calculation process, and can effectively improve the calculation efficiency via a segmented calculation in the time dimension. The calculation results are in good agreement with the traditional analytical solution, indicating that the method is reasonable and feasible. b. The method in this paper considers the case in which the external load is a function of time and can solve the calculation of the consolidation degree and the settlement of the foundation under the conditions of an impulse load, a simple harmonic load, or multiple step loading in actual engineering. Only the function expression of the external load is input into the program during the calculation, which avoids the steps of other theoretical calculation methods that need to convert the external load. Via this method, human error is reduced, and the calculation process for the degree of consolidation and the settlement of the foundation is greatly simplified, with good engineering application value. c. The application condition of the presented method is the one-dimensional consolidation of a single-layered saturated soil with a permeable surface. The calculation formulas of the finite difference method should be transformed when analyzing the one-dimensional consolidation of the multilayered unsaturated soil under partially permeable boundary conditions. Besides, a multidimensional consolidation issue cannot be solved by this method.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.