The Terzaghi theory for one-dimensional consolidation states that all quantifiable changes in the stress of a soil (compression, deformation, and shear resistance) are a direct result of changes in effective stress. The effective stress σ ′ is related to total stress σ and the pore pressure p by the following relationship: σ = σ ′ + p (1)

The overpressure dissipation is described by the following equation: C v ∂ 2 u ∂ z 2 = ∂ u ∂ t , (2) C v = k E s γ w , (3) where C v is the coefficient of consolidation, E s is the modulus of compressibility, k is the hydraulic conductivity, γ w is the unit weight of water, and u is the excess pore pressure.

In addition, Terzaghi’s consolidation theory (Terzaghi, 1929; Terzaghi, 1943) was put forward under the following assumptions:

(1) The soil is homogeneous.

It is therefore necessary to point out that according to the principle of effective stress, total stress increment is produced by the load applied on the multilayer systems. Resultantly, the effective stress increases with excess pore pressure decrease as time goes on. σ ′ = σ − P (4)

According to Eq. 4, the following equation with the variable of effective stress increment can be obtained, which will benefit our solution. C v ∂ 2 σ ′ ∂ z 2 = ∂ σ ′ ∂ t (5)

The common weighted average method generates a C v for the whole multilayer system according to C v = ∑ i = 1 n h i × C v i ∑ i = 1 n h i . (6)

The Laplace transform is a widely used integral transform with many applications in physics and engineering. Hence, Stehfest numerical inversion of Laplace transforms is adopted (Stehfest, 1960) for our study, and we utilize the transformation matrix to connect parts within the multilayer system. According to the integral multilayer transformation matrix and transformation matrix between different layers, the distribution of effective stress increments and excess pore pressure can be obtained. Moreover, results are verified with an implicit finite difference numerical solution. Figure 2 shows a schematic plot of the multilayer Terzaghi consolidation.

Consequently, the Laplace transform of Eq. 5 can be expanded to C v ∂ 2 σ ˜ ′ ( z , s ) ∂ z 2 = s σ ˜ ′ ( z , s ) −   σ ′ ( z , 0 ) , (7) where σ ˜ ′ ( z , s ) is the Laplace transform of σ ′ ( z , t ) .

At the beginning of consolidation, according to the effective stress principle, pore pressure is equal to overburden stress, which is zero initial effective stress; thus Eq. 7 can be obtained. C v ∂ 2 σ ˜ ′ ( z , s ) ∂ z 2 = s σ ˜ ′ ( z , s ) . (7a)

Then, the general solution of the ordinary differential equation for Eq. 7 is as follows: σ ˜ ′ ( z , s ) =   c 1 ⁡ exp ( β z ) + c 2 ⁡ exp ( − β z ) , (8) where c 1 and c 2 are constants, β = s c v .

Combining Eq 9 错误!未找到引用源。 and its partial derivative about z, the following expression can be derived: [ σ ′ ∼ ( z , s ) ∂ σ ′ ∼ ( z , s ) ∂ z ] = [ exp ( β z ) exp ( − β z ) β ⁡ exp ( β z ) − β ⁡ exp ( − β z ) ] [ c 1 c 2 ] . (9) When z = 0 : [ σ ′ ∼ ( 0 , s ) ∂ σ ′ ∼ ( 0 , s ) ∂ z ] = [ 1 1 β − β ] [ c 1 c 2 ] . (9a)

Combining Eq 10 and Eq 11 yields Eq. 12: [ σ ′ ∼ ( 0 , s ) ∂ σ ′ ∼ ( 0 , s ) ∂ z ] = [ 1 2 [ exp ( − β z ) + exp ( β z ) ] 1 2 β [ exp ( − β z ) − exp ( β z ) ] β 2 [ exp ( − β z ) − exp ( β z ) ] 1 2 [ exp ( − β z ) + exp ( β z ) ] ] [ σ ∼ ′ ( z , s ) ∂ σ ′ ∼ ( z , s ) ∂ z ] . (10) And when z i is not zero: [ σ ′ ∼ ( z i , s ) ∂ σ ′ ∼ ( z i , s ) ∂ z ] = [ 1 1 β − β ] [ exp ( β z i ) 0 0 exp ( − β z i ) ] [ β 2 β ⁡ exp ( β z ) 1 2 β ⁡ exp ( β z ) β 2 β ⁡ exp ( − β z ) − 1 2 β ⁡ exp ( − β z ) ] [ σ ′ ∼ ( z , s ) ∂ σ ′ ∼ ( z , s ) ∂ z ] . (11)

Then the relationship between the top surface stress and the bottom stress can be derived. When considering the equation of continuous stress, and flow conservation, between two layers, the relationship between different layers can be derived. k i ∂ σ ′ ∼ ( z i − , s ) ∂ z = k i + 1 ∂ σ ′ ∼ ( z i + , s ) ∂ z , (12) σ ′ ∼ ( z i − , s ) = σ ′ ∼ ( z i + , s ) . (13)

Combining Eq 12 and Eq 13 yields Eq. 16: [ σ ′ ∼ ( z i − , s ) ∂ σ ′ ∼ ( z i − , s ) ∂ z ] = [ 1 k i + 1 k i ] [ σ ′ ∼ ( z i + , s ) ∂ σ ′ ∼ ( z i + , s ) ∂ z ] . (14)

Stress distribution in the same layer and in the interface can be derived, respectively, by using Eq. 11 and Eq. 14. With the equation of each layer combined together, a transform matrix T can be obtained to express the relationship between z = 0 and z = z n . [ σ ′ ∼ ( 0 , s ) ∂ σ ′ ∼ ( 0 , s ) ∂ z ] = T [ σ ′ ∼ ( z n , s ) ∂ σ ′ ∼ ( z n , s ) ∂ z ] (15)

Here, this article only considers the situation of the pervious top surface and impervious bottom PTIB for illustration.

The following boundary conditions are implemented: z = 0,u ( z,t ) = 0;  z = z n  ,  ∂ u ( z , t ) ∂ z = 0. (16)

The corresponding Laplace transformation: z = 0 , σ ′ ∼ ( 0 , s ) = σ s ; z = z n , ∂ σ ′ ∼ ( z n , s ) ∂ z = 0 (17) Hence σ ′ ∼ ( z n , s ) = σ T 11 S , (18) where σ is the pressure on the surface and T₁₁ is the value of the first column and the first row of T.

With σ ′ ∼ ( z n , s ) , the stress at each upper point σ ′ ∼ ( z , s ) can be obtained by the transformation matrix. Moreover, the real stress distribution can be derived by the inverse of Laplace transformation.

As for an n-layer consolidation system, the multilayer consolidation transform matrix is T 1 = [ 1 2 [ exp ( − β 1 h 1 ) + exp ( β 1 h 1 ) ] k 2 k 1 1 2 β 1 [ exp ( − β 1 h 1 ) − exp ( β 1 h 1 ) ] β 1 2 [ exp ( − β 1 h 1 ) − exp ( β 1 h 1 ) ] k 2 k 1 1 2 [ exp ( − β 1 h 1 ) + exp ( β 1 h 1 ) ] ] T i = [ 1 2 [ exp ( − β i h i ) + exp ( β i h i ) ] k i + 1 k i 1 2 β i [ exp ( − β i h i ) − exp ( β i h i ) ] β i 2 [ exp ( − β i h i ) − exp ( β i h i ) ] k i + 1 k i 1 2 [ exp ( − β i h i ) + exp ( β i h i ) ] ] T n = [ 1 2 [ exp ( − β n h n ) + exp ( β n h n ) ] 1 2 β n [ exp ( − β n h n ) − exp ( β n h n ) ] β n 2 [ exp ( − β n h n ) − exp ( β n h n ) ] 1 2 [ exp ( − β n h n ) + exp ( β n h n ) ] ] [ σ ′ ∼ ( 0 , s ) ∂ σ ′ ∼ ( 0 , s ) ∂ z ] = T 1 ⋅ ⋅ T i ⋅ ⋅ T n . [ σ ′ ∼ ( z n , s ) ∂ σ ′ ∼ ( z n , s ) ∂ z ] ,   i = 2 , 3 ....... n − 1 (19)

The common weighted average method will lead to the following weighted average method transform matrix: [ σ ′ ∼ ( 0 , s ) ∂ σ ′ ∼ ( 0 , s ) ∂ z ] = [ 1 2 [ exp ( − β z n ) + exp ( β z n ) ] 1 2 β [ exp ( − β z n ) − exp ( β z n ) ] β 2 [ exp ( − β z n ) − exp ( β z n ) ] 1 2 [ exp ( − β z n ) + exp ( β z n ) ] ] [ σ ′ ∼ ( z n , s ) ∂ σ ′ ∼ ( z n , s ) ∂ z ] . (20)

When t → ∞ , s → 0 , β → 0 , by applying the following Taylor expansion we obtain e x = ∑ n = 0 ∞ x n n ! . (21)

As for a multilayer transform matrix, Eq. 18 will develop into [ σ ′ ∼ ( 0 , s ) ∂ σ ′ ∼ ( 0 , s ) ∂ z ] = [ 1 − ( h n + ∑ i = 1 n − 1 k n k i h i ) 0 k n k 1 ] [ σ ′ ∼ ( z n , s ) ∂ σ ′ ∼ ( z n , s ) ∂ z ] . (22)

Then, the weighted average method transform matrix given in Eq. 20 will develop into [ σ ′ ∼ ( 0 , s ) ∂ σ ′ ∼ ( 0 , s ) ∂ z ] = [ 1 − z 0 1 ] [ σ ′ ∼ ( z , s ) ∂ σ ′ ∼ ( z , s ) ∂ z ] . (23)

It can therefore be concluded that when conductivities in different layers are nearly the same, the weighted average method can be adopted for the whole multilayer system, which is the situation of homogenous.

This analytical solution is verified against the results of Lee et al. (1992), as shown in Figure 3. The parameters used in Lee’s model are shown in Table 1, and the model sketch map follows the illustration in Figure 1 (right), showing a simplified sketch map of four layers with a PTIB and without overlaying water.

To compare different upscaling methods, an implicit finite-difference numerical model is developed. This study utilizes a numerical model as a benchmark model and is verifiable with both laboratory tests and a basin model (Jingchen Zhang., 2015). Analytical results are also compared with implicit finite difference numerical solution for verification. The parameters of the three different layers are shown in Table 2. Comparison of the results is then shown in Figure 4. Hence, a conclusion can be reached that this method can be applied to multilayer Terzaghi consolidation compaction. It should be noted that the values of m v l and surcharge ensure the small strain, and correspondingly nearly no settlement can be seen in Figure 4. Resultantly, the small changes in height are shown in Figure 5.

Here, T represents the properties of the multilayer system; then a new consolidation coefficient C v is needed to represent the whole multilayer system. With the new C v , a new transform matrix for the multilayer consolidation T ′ can be obtained. The new C v should be the one that makes the minimum of Eq. 24, to numerically match T and T ′ . ( T ′ ( 1,1 ) − T ( 1,1 ) ) 2 + ( T ′ ( 1,2 ) − T ( 1,2 ) ) 2 + ( T ′ ( 2,1 ) − T ( 2,1 ) ) 2 + ( T ′ ( 2,2 ) − T ( 2,2 ) ) 2 (24)

Under the condition of PTIB, it can be found that ∂ σ ′ ∼ ( z n , s ) ∂ z = 0 , so σ ∼ ( 0 , s ) = T ( 1,1 ) × σ ′ ∼ ( z n , s ) ′ . This provides a thought of using a homogeneous layer T ( 1,1 ) to represent multilayer nonhomogeneous consolidation. A new β is required to fit the value of T ( 1,1 ) . The example of the afore-mentioned three layers with different properties is adopted for this illustration. The new C v can be derived from 1 2 [ e − β z + e β z ] = T ( 1,1 ) . (25)

C v changes with time according to Eq. 25, as shown in Figure 6.

In light of the long-time consolidation, C v is set to be 1.4 × 10⁻⁸ (m²/s). A homogeneous layer with the new upscaling C v can then be compared with the three-layer system. Specifically, the transform matrix for a three-layer system is T 1 = [ 1 2 [ e − β 1 h 1 + e β 1 h 1 ] k 2 k 1 1 2 β 1 [ e − β 1 h 1 − e β 1 h 1 ] β 1 2 [ e − β 1 h 1 − e β 1 h 1 ] k 2 k 1 1 2 [ e − β 1 h 1 + e β 1 h 1 ] ] T 2 = [ 1 2 [ e − β 2 h 2 + e β 2 h 2 ] k 3 k 2 1 2 β 2 [ e − β 2 h 2 − e β 2 h 2 ] β 2 2 [ e − β 2 h 2 − e β 2 h 2 ] k 3 k 2 1 2 [ e − β 2 h 2 + e β 2 h 2 ] ] . T 3 = [ 1 2 [ e − β 3 h 3 + e β 3 h 3 ] 1 2 β 3 [ e − β 3 h 3 − e β 3 h 3 ] β 3 2 [ e − β 3 h 3 − e β 3 h 3 ] 1 2 [ e − β 3 h 3 + e β 3 h 3 ] ] [ σ ′ ∼ ( 0 , s ) ∂ σ ′ ∼ ( 0 , s ) ∂ z ] = T 1 T 2 T 3 . [ σ ′ ∼ ( z 3 , s ) ∂ σ ′ ∼ ( z 3 , s ) ∂ z ] . (26)

Meanwhile, the common weighted average method will lead to another C v . Here, C _v = 3.0093 × 10⁻⁸ (m²/s) in the weighted average method is used. With this C v , a new β and a weighted average transform matrix can be derived. Comparisons between the three methods are shown in Figure 7. R = ∑ i = 1 n | u i ′ − u i | / u i n (27)

According to Eq. 27, the relative error between the weighted average method or the new upscaling method and the three-layer numerical results can be obtained, as shown in Figure 8.

When it comes to 20,000 days, there is nearly no overpressure; hence, the multilayer system’s characteristics are studied within 20,000 days. As can be seen from Figure 9, in the first 100 days, the weighted average method is more efficient than the upscaling method. As fluid flows out through the top surface, the whole system is determined by the properties of the first layer before pressure reduction reaches the second layer. Overpressure distributions of one homogeneous layer (same properties with the upper layer) and three nonhomogeneous layers’ consolidation are the same. This can explain why the results of changing C v and the weighted average method compact faster than the real situation. The possible explanation is that bigger C v of the upper layer is applied to the whole layers with small C v characteristics. The multilayer consolidation will show the integral properties more accurately after the stimulation reaches the bottom.

Gibson et al. (1967) developed a large-strain consolidation theory with more general basic assumptions than the small-strain theory in (3), (5), and (7) of Terzaghi’s assumptions. The limitation of small strains was not imposed, and the soil compressibility and permeability are allowed to vary with the void ratio during consolidation. These assumptions are closer to the actual scenario. Furthermore, Darcy’s law is assumed to be valid, but it is recasted in a form in which it is the relative velocity of the soil skeleton and the pore fluid. The fluid velocity is related to the excess pore pressure gradient. − 1 1 + e 0 ∂ e ∂ t + ( γ s γ w − 1 ) d d e [ k 1 + e ] ∂ e ∂ z + ∂ ∂ z [ − k ( 1 + e 0 ) γ w ( 1 + e ) d σ ′ d e ∂ e ∂ z ] = 0 , (28) where e is the void ratio; e 0 is the initial void ratio; k is the conductivity; γ w and γ s are the unit weight of water and soil, respectively; z is the solid coordinate; and σ ′ is the effective stress.

To apply the transfer matrix method to large-strain consolidation, the simplification of Xie et al. is adopted to simplify Gibson’s equation (Xie and Leo, 2004). The influence and theoretical analysis of this simplification can be found in his research.

The coefficient of volume compressibility of the soil skeleton for a large strain is defined as − 1 1 + e d e d σ ′ = m v l . (29)

Then Gibson’s equation can be changed to 1 γ w ∂ ∂ z [ k ( 1 + e 0 ) ( 1 + e ) ∂ u ∂ z ] = m v l 1 + e 1 + e 0 ( ∂ u ∂ t − ∂ q u ∂ t ) , (30) where q u is the surcharge and m v l is constant during consolidation.

The relationship between conductivity k and void ratio is k k o = ( 1 + e 1 + e 0 ) 2 , (31) where k o is the initial conductivity at time t = 0, k is often found empirically to be a logarithmic function of the void ratio, and e₀ is the initial void ratio.

A load q u is applied suddenly at t = 0 on the top surface of the model and remains constant thereafter. According to the effective stress principle and Xie’s assumption (Xie, K. & Leo, C., 2004), the relationship between the void ratio and the excess pore water pressure can be deduced as follows. 1 + e 1 + e 0 = exp ( − m v l ( q u − u ) ) (32)

With Eq. 29, Eq. 30 and Eq. 32 can now be changed to the following one, which determines excess pore evolution in Lagrangian coordinates: c v o [ ∂ 2 u ∂ z 2 + m v l ( ∂ u ∂ z ) 2 ] = ∂ u ∂ t , (33) where c v o is the initial coefficient of consolidation at time t = 0 given by c v o = k o m v l r w . (34)

The solution to the large consolidation theory is facilitated by the following transformation. ω = ω ( z , t ) = exp ( m v l u ) (35)

In consideration of the permeable top impermeable base (PTIB) boundary condition, with Eq. 35 and Eq. 33, the equation of the Terzaghi form will be obtained as follows: c v ∂ 2 ω ∂ z 2 = ∂ ω ∂ t ω ( 0 , t ) = 1 ∂ ω ∂ z ( H , t ) = 0 ω ( H , 0 ) = exp ( m v l q u ) . (36)

Then the transfer matrix can be adopted for the solution and upscaling of multilayer Gibson consolidation. To use the transfer matrix method, ω should be continuous; hence, different layers share the same m v l . Then according to Eq. 34, what upscaled is actually k o , i.e., conductivity upscaling.

Analytical solution is compared with the developed implicit finite difference numerical code for verification. The model sketch map follows Figure 1 错误!未找到引用源。 (right), PTIB, and no overlying water. A three-layer model is adopted for comparison, with the surcharge 1 × 105 Pa, γ _s = 26.950 kN/m³, and γ _w = 10.045 kN/m³. For each layer, m _vl = 4 × 10⁻⁶ Pa⁻¹ and thickness = 5 m. The k o of the surface, middle, and bottom layers are 1.00 × 10⁻⁹ m/s, 1.16 × 10⁻¹⁰ m/s, and 1.04 × 10⁻⁹ m/s. The e 0 of the surface, middle, and bottom layers are 3, 4, and 5. Figure 10 shows the comparison results, which are consistent and proves the effectiveness of this method in solving large-strain multilayer consolidation.

To evaluate this upscaling method, comparison with the weighted average method is carried out here, utilizing the same three-layer model in Verification. When it comes to 1,000 years, there are nearly no excess pore pressures according to the weighted average method; hence, we focus within 1,000 years.

Upscaling k v o changes with time as shown in Figure 11. In consideration of the long geology process, the k o value is set to be 2.8356 × 10⁻¹⁰ (m/s), and C v is 7.2337 × 10⁻⁹ (m²/s). While for the weighted average method, k o is 7.1867 × 10⁻¹⁰ (m/s) and C v is 1.8333 × 10⁻⁹ (m²/s). The comparisons of the two upscaling methods and numerical solution are shown in Figure 12.

The whole model only shows properties of the first layer before the pressure reduction reaches the second layer. Also, the weighted average k o is closer to the first layer’s k o than the upscaled k o . Hence, within the first 10 years, the results obtained through the weighted average method are closer to the numerical results than those obtained through the upscaling method. However, as a whole, the upscaling method is more effective than the weighted average method.

The multilayer system only shows properties of the place affected by the stimulation, and the integral properties will change with the increase in the affected region. This can partly explain the result of a changing upscaled k o . Consequently, this upscaling method is more accurate than the common weighted average method as a whole for multilayer large-strain consolidation.