Nonlinear consolidation analysis of multilayered soil with coupled vertical-radial drainage using the spectral method

The nonlinear variation of soil compressibility and permeability with void ratio (i.e., e-log σ′ and e-log k) has been included in the consolidation theory to accurately predict the behavior of soft soil stabilized by vertical drains. However, most current nonlinear consolidation models incorporating the coupled radial-vertical flow are based on some simplified assumptions, while including some features such as the complex implementation of multilayered computations, time-dependent loading and stress distribution with depth. This study hence introduces a novel approach where the spectral method is used to analyze the nonlinear consolidation behavior of multilayered soil associated with coupled vertical-radial drainage. In addition, time- and depth-dependent stress and soil properties at each soil layer are incorporated into the proposed model. Subsequently, the solution is verified against experimental and field data with comparison to previous analytical solutions. The results show greater accuracy of the proposed method in predicting in-situ soil behavior. A parametric study based on the proposed solution indicates that the ratio between the compression and permeability indices (ω = Cc/Ck) has a great impact on the consolidation rate, i.e., the greater the ω, the smaller the consolidation rate. Increasing the load increment ratio and the absolute difference between unity and ω (i.e., |ω − 1|) can exacerbate prediction error if the conventional simplified methods are used.

Note that the above features can be combined to provide improved predictions. However, most of them were based on simplified assumptions of constant soil compressibility and permeability during consolidation.
As the sedimentary history and stress conditions of soil can vary significantly in the field, most soft soils are rarely homogeneous and usually consist of several layers [2]. However, previous nonlinear consolidation models show limited capacity in capturing the influence of adjacent soil layers because they strictly rely on specific loading and stress distribution patterns. This study, therefore, aims to overcome the above limitations in previous studies [71,72,74,78] by considering the nonlinear compressibility and permeability based on the spectral method framework, so that a more realistic and rigorous solution for the PVD-assisted soil consolidation can be achieved. In this paper, the spectral method is adopted to solve the governing equations, and subsequently, the model is verified against the experimental and field data in comparison with previous simplified solutions. Finally, the applicability and threshold limits of the past and the current solutions are discussed and evaluated.

Limitations of existing models
This section firstly details the limitations of existing mathematical solutions, followed by the objectives and innovations of the current study. The logarithmic models (e'-log r and e-log k) are commonly used to represent the variations of soil compressibility and permeability with void ratio, which can be represented by [64]: where e is the void ratio while the subscript 0 denotes the initial state; C c , C r , C kh , C kv are the compression index, the recompression index, the radial permeability index and the vertical permeability index, respectively; r 0 0 , r 0 p and r 0 are the initial effective stress, the yield stress (effective preconsolidation pressure) and the average effective stress, respectively; k h and k v are the radial and vertical permeability coefficients of the undisturbed soil, respectively.
Overall consolidation degree The following parameters are now introduced and defined as: Then the radial and vertical consolidation coefficients can be expressed as: The above expressions (i.e., Eqs. (4)-(8)) show how the compression and permeability of soil would change due to the reduced void ratio during consolidation. Due to the complexity in solving the consolidation governing equations, some studies [16,17,30] assumed that B = B = -1 based on the field situation, where the compression C c is very close to the permeability indices (C kv or C kh ), while the others (summarized in Table 1) have obtained simplified analytical solutions based on the following assumptions: (1) Simplified Method A: use an average value to represent the ratio of the effective stress to the initial effective stress (i.e., r 0 r 0 0 in Eqs. (7) and (8)) during the consolidation process, i.e., r 0 where q(t), q max and u are the time-dependent loading, the final level of loading and EPWP, respectively [46,47,75]; (2) Simplified Method B: use the average values to represent the varying consolidation coefficients, which are the nonlinear coefficient terms in the governing equation, i.e., s h i [28,32,33,35,38,41,79]. Table 1 lists the capabilities and assumptions of some significant nonlinear consolidation models. It can be seen from Table 1 that the main limitations of previous nonlinear consolidation models are as follows: (a) Although Simplified Methods A and B based on Assumptions (1) and (2) adopt the voidratio-stress relationship (Eq. (1)) for settlement and EPWP calculations, these two assumptions make the consolidation coefficients (i.e., the coefficient terms of the consolidation governing equation) constant. This means the nonlinear behavior is not included in the dissipation equation of EPWP properly [28,32,33,35,38,41,46,47,75].
Þln 10 for r 0 p \r 0 (b) Simplified Methods A and B directly adopt these assumptions to linearize the nonlinear coefficient terms of the governing equation. The validity and associated threshold have not been established. In other words, the acceptance range of error caused by the simplified assumptions has not been evaluated. It is necessary to evaluate the errors induced by simplified assumptions that help understand the validity of Simplified Methods A and B, and thus determine the appropriate range of soil parameters [28,32,33,35,38,41,46,47,75]. (c) Although some of the nonlinear consolidation models can consider coupled radial-vertical drainage, they can only consider a single layer of soil while changes in soil parameters and stress distribution along the depth are neglected [41,47,75].
In view of the above, the objectives of this study are to provide a more general nonlinear consolidation model which can consider the following factors: i. Coupled vertical-radial drainage; ii. Nonlinear permeability and compressibility during consolidation process; iii. Multilayered condition; iv. Time-dependent loading; v. Over-consolidated and normally consolidated state.
The key advantages of the current approaches are shown in Fig. 2 in comparison with conventional methods.

Basic assumptions
The following basic assumptions in this study were adopted while developing the mathematic model.
(a) The soil particles and water are incompressible. The nonlinear relationships between void ratio with permeability and effective stress during consolidation are shown in Eqs. (2) and (3). (b) The compressibility and the vertical permeability coefficients in the smear and undisturbed zones are assumed to be the same. The horizontal permeability coefficient in the smear zone is constant distribution and the ratio of horizontal permeability coefficients outside and in the smear zone is constant during consolidation. The size of the smear zone is constant throughout the depth. (c) The initial effective stress, the pre-consolidation pressure, the vertical stress, and associated parameters for a given lth layer of soil with relatively small thickness are assumed to be constant, but they change with depth as shown in Fig. 3. (d) The soil is assumed to be fully saturated, and the velocity of pore water flow is governed by Darcy's law. Although the EPWP varies in the radial direction, the average EPWP along the radial direction is used at a given depth to combined with flow in the vertical direction as shown in Eq. (9), following the approach proposed by Tang and Onitsuka [62]. (e) Strains only occur in the vertical direction, which are equal at a given depth along the radial direction (equal strain condition).

Governing differential equations
The unit cell for a multilayered soil with a vertical drain is shown in Fig. 3. The governing differential equation for soil consolidation, while considering vertical and radial drainage, can be given by (see ''Appendix A1'' for derivation): where u is the average EPWP at a particular depth; t is the time; H is the total depth of soil; Z is the normalized depth, i.e., Z¼z=H in which z is the depth; c w is the unit weight of water; m v0 is the initial volume compressibility, and it can be calculated by when r 0 r 0 p ; l is the dimensionless parameter, which is computed based on the permeability variation of soil within the smear zone, the radial geometry of the drain. Detailed calculation of l can be referred to the previous studies, e.g., Walker and Indraratna [71], Lu et al. [47] and Nguyen [49].
Given a time-dependent loading q t ð Þ, the effective stress can be determined by: the governing equation can be rewritten as: It can be seen from Eq. (11) that when vertical drainage is not considered and q t ð Þ becomes the instantaneous loading, the above equation turns into the nonlinear radial consolidation model of Walker et al. [73] without considering non-Darcian flow. If the above nonlinear term is further replaced by the average value (i.e., , it becomes the nonlinear radial consolidation model by Indraratna et al. [28]. Furthermore, when B h-= 1 (the slopes of e-log r 0 and e-log k are the same, the above governing equation becomes the same as that by Hansbo [20].

Advanced features of spectral method
Spectral method is one of the very advanced mathematical techniques for facilitating numerical solution of even complex partial differential equations (PDEs). It evolved after the common numerical category of finite element method (FEM) and finite difference method (FDM) whose the accuracy depends on the size of the subdomain. [67]. The spectral method is based on global basis functions (highorder polynomial or trigonometric functions). Compared with the numerical methods such as FDM and FEM, the spectral method has the following advantages when the geometry of the problem is fairly smooth and regular (e.g., consolidation) [6,52]: (1) high calculation accuracy; (2) memory-minimizing and computational efficiency; and (3) high stability. Therefore, the method can capture the transition of variables over time and space such as stress, EPWP and soil properties. It was adopted in the current study to solve the complex governing equation incorporating the variation of multiple soil properties during consolidation. When the pore pressure profile changes sharply, oscillations may occur near steep fronts, which is called the Gibbs phenomenon. The Gibbs phenomenon can be reduced or eliminated by increasing the series of N term. Therefore, more series terms are required when modeling sharp changes in the pore pressure profile.

Solutions based on the spectral method
For the spectral method, the EPWP u Z; t ð Þ is expressed as a truncated series of N terms, which can be expressed in matrix form as follows [71,72,74,78]: where U j are known basis functions and A j are expansion coefficients which can vary with time, and U¼ U 1 U 2 ::: U N ½ ð 13Þ The choice of the basis functions needs to satisfy the boundary conditions of governing equation [6]. The pervious top-pervious bottom (PTPB) and the pervious topimpervious bottom (PTIB) boundary conditions are given, respectively, by: With respect to these boundary conditions (Eqs. (15) and (16)), the appropriate choice of basis function can be given by: Note that in the current study, the material properties such as permeability and compressibility vary with void ratio and the effective stress, resulting in the complexity to obtain an exact solution through the spectral method. Therefore, the current study proposes a numerical approach where the consolidation process is divided into a discrete number of time steps (Fig. 4). During each time step, the material parameters are assumed to be constant, but they are then re-computed and updated in the next time step based on Eqs. (6)- (8). By updating the material properties at each time step and combining the weighted residual method (WRM), A t ð Þ can be obtained using Eq. (19), thereby the EPWP at a given depth and time can be obtained in a matrix form (see ''Appendix A2'' for derivation): where the elements of matrices C, I and W incorporate the loading patterns and material parameters of every soil layer. The detailed expressions of these elements can be found in ''Appendix A2'', Eqs. (40)- (43). Figure 4 is the flowchart showing the detailed implementation of the proposed model. Note that the interval of time step affects the accuracy of the updated soil parameters for the next time step.
Since the EPWP at a given depth is expressed as a function of depth and time as shown in Eq. (20), the average pore water pressure u avg Z l ; Z lþ1 ; t ð Þand settlement S Z l ; Z lþ1 ; t ð Þin the lth layer (between depths Z l and Z l?1 ) can be calculated, respectively, by: The overall average degrees of consolidation for the multilayered soil defined by the excess pore pressure and settlement can, respectively, be obtained by: where Z l and Z l?1 denote the normalized depth at the bottom and top of the lth layer, respectively; K l s is the stress influence factor in the lth layer; S Z l ; Z lþ1 ; 1 ð Þis the final settlement in the lth layer. The superscript l represents the value of the corresponding parameter or variable at the lth layer.

Continuity conditions at the soil interface
In the solution of the spectral method, the average EPWP u Z; t ð Þ is expressed as a truncated series of N terms, and the sine functions were selected as the basis functions. Therefore, the value of the average EPWP and its derivative in the soil at any position are continuous.
First, the continuous condition of EPWP at the interface between two adjacent layers (lth and l ? 1th layer) can be satisfied: In addition, since the soil property of a certain soil layer is assumed to be constant in this study, an interface (i.e., dummy) layer with a thickness of zero is set between two adjacent layers, as shown in Fig. 3. The soil parameters are assumed to be linearly distributed in the interface layer, so the continuous condition of flow rate can also be satisfied between two adjacent layers (lth and l ? 1th layer), i.e., It is noteworthy that the distribution made by the dummy layer to C ij and I i is zero, and the distribution made by the interface layer between two adjacent layers (l th and l ? 1 th layer) to W ij can be found in Eq. 43.

Model verification
To verify this proposed model, the mathematical formulation presented above is applied to the following laboratory and field studies: 1. Radial consolidation of single soil layer by singleinstantaneous loading [28,73]; 2. Vertical and radial consolidation of multilayered soil by multi-ramp loading [10,60]; The calculation of the dimensionless drain parameter l is based on the assumption of a smear zone with constant reduced permeability [20].

Laboratory tests
Two laboratory studies [28,73] were used to verify the proposed model. The physical size of the consolidation apparatus was 450 mm in diameter by 950 mm high, and the reconstituted alluvial clay from Moruya (New South Wales) was used. For these tests (normal consolidation range), the initial pre-consolidation pressures r 0 0 of the soil were 20 kPa and 50 kPa with the loading increments in these two studies were 30 kPa and 50 kPa, respectively (i.e., Dp = 30, 50 kPa). The detailed testing procedure can be found in Walker et al. [73]. The soil parameters and drain properties are shown in Table 2. Note that as the drain was relatively short, the well resistance effect was neglected in the calculation of l.
The degree of consolidation based on the settlement is obtained using Eq. (24). The accuracy of the calculation is determined by the selection of the truncated series N, as shown in Eq. (12). An investigation on the convergence was carried out especially addressing the effects of the numbers of the truncated series N through these two Pre-consolidation pressure, r 0 0 (r 0 p )/kPa 20 50 Load, p/kPa 30 50 laboratory tests, the results are shown in Fig. 5. It shows that N = 50 are sufficient for a single soil layer with an error \ 0.5% for calculating EPWP. In addition, Fig. 5b shows the ''exponential'' convergence of spectral method with N. The relationship between N and error d can be expressed as log(N) = a ? b log(d), where a and b are the coefficients. In practical applications, the selection of N depends on the complexity of the problem (e.g., the number of soil layers and the differences in parameters between soil layers) and the required accuracy of computation. When the number of soil layers is large or the soil parameters differ greatly, the value of N should be larger to improve the calculation accuracy and eliminate the influence of the Gibbs phenomenon [6,52]. The appropriate truncation series N can be selected according to the calculation accuracy requirements. For example, if the number of digit accuracy (p) is required, the truncated series N can be selected based on the relationship of N ¼ 10 aþb logð10 Àp Þ in this case. In these two tests, only radial drainage was allowed, so dT v0 was set as 0. The results were compared with the laboratory data and the analytical solutions presented by Indraratna et al. [28], Walker et al. [73] and Lu et al. [46], as shown in Fig. 6a. Note that the model of Walker et al. [73] is a closed-form analytical solution, while the models of Lu et al. [46] and Indraratna et al. [28] are simplified analytical solutions based on Simplified Methods A and B, respectively. Figure 6 shows that the results calculated by the proposed method are very close to the experimental results and analytical solution of Walker et al. [73]. Indeed, Fig. 6b and c shows that the largest deviations between the proposed solution and measured data and analytical solutions [73] in Test 1 and Test 2 are less than 4.6% and 0.6%, respectively. The difference between the calculation results of the proposed method and the analytical solution of Walker et al. [73] is caused by the insufficient value of truncation series N. When the value of N increases, the result predicted by the proposed method becomes more accurate and closer to the analytical solution. However, the results from the models of Lu et al. [46] and Indraratna et al. [28] deviate from accuracy, especially in the early stage of consolidation. This is because the average consolidation coefficients have been in these two models, which overestimate the actual consolidation coefficient during the early stage, as shown in Fig. 6d and e.

Hangzhou-Ningbo (HN) Expressway, China
The test embankment using PVDs at Hangzhou-Ningbo (HN) Expressway was reported by Chai et al. [10] and Shen et al. [60]. The HN Expressway was located at the southern coast of Hangzhou Bay, China. The thickness of the soft layers was about 23 m. The top crust was considered to be in lightly over-consolidated state with an over-consolidation ratio (OCR) of about 5, and the deeper layers were in the normally consolidated state. The soil profile and soil parameters used in this study provided by Chai et al. [10] are shown in Table 3. The stage loading process is shown in Fig. 7a. The final fill height for the surcharge preloading was 5.88 m and the unit weight of the fill material was 20 kN/m 3 . As suggested by Tavenas et al. [64], the permeability indices in this study were calculated by C kh = C hv = 0.  [32,46,74]. Note that the model of Walker et al. [74] is the conventional linear consolidation model for multilayered  soil with coupled vertical-radial drainage, the initial permeability coefficients (k h0 and k v0 ) and compressibility coefficient (m v ) were adopted in the linear consolidation model, as shown in Table 3. The settlements predicted by the model of Walker et al. [74] overestimate the field data due to the inability to consider the nonlinear behavior of the soil. The models of Lu et al. [46] and Kim et al. [32] are analytical solutions of radial consolidation (i.e., only radial drainage) based on the Simplified Methods A and B, respectively, which can be considered as the piecewise solutions in which the soil parameters and stress conditions are the corresponding values at the mid-point of each layer.
Since the ratios between the compression and permeability indices of the main compression soil layer (e.g., 4.8-19 m) are close to 1 in this case (see Table 3), the settlements predicted by the models of Lu et al. [46] and Kim et al. [32] are very similar. However, since the nonlinear vertical permeability are not included in these two models, their predicted settlements underestimate the field data. The proposed method which incorporates both nonlinear vertical and radial permeability provides better prediction. For example, the difference between the predicted and measured settlement at 400 days significantly decreases to about 10 mm (0.58%) using the proposed model, where the errors in the analyses of Walker et al. [74], Lu et al. [46] and Kim et al. [32] are 71 mm (4.11%), 82 mm (4.73%) and 87 mm (5.02%), respectively. Figure 8 compares the predicted EPWPs by the current proposed solution with the measured data and the results obtained by the methods of Walker et al. [74], Lu et al. [46] and Kim et al. [32] at three different depths (i.e., z = -2.0 m, z = -10.0 m and z = -14.05 m). The EPWPs predicted by the model of Walker et al. [74] overestimate the dissipation rate of EPWPs at all depths, as this model cannot consider the nonlinear behaviors of the soil. Generally, the results by the proposed method are closer to the field data compared to other models, especially in shallow soil, i.e., at depth of 2 m. For example, at 200 days, the error of 30 kPa in previous models is reduced to be less than 3 kPa by the proposed model. At greater depths (i.e., z = -10.0 m and z = -14.05 m), the ratios between the compression and permeability indices (C c /C kh and C c /C kv ) are close to 1, and the soil consolidation is predominantly governed by the radial drainage. Therefore, the EPWPs predicted by models of Lu et al. [46] and Kim et al. [32] approach closer to the field data and the current model, as shown in Fig. 8b and c. Note that all the predicted EPWPs dissipate completely after 800 days while the measured EPWPs gradually change after 600 days. For example, the measured EPWP at 2 m depth remains almost unchanged at about 10 kPa until the end of observation. This residual EPWP could be attributed to the effect of rising groundwater level after 600 days. Figure 9 shows the distribution of EPWP along the depth at 100, 200 and 400 days. In general, the isochrones of EPWP are in good agreement with the measured data. In fact, the predicted curves present very well the smooth transition in EPWP over different soil layers, provided appropriate value of N. This shows that the proposed model based on the spectral method can be well applied to the nonlinear consolidation calculation of multilayered foundations.  Fig. 10 Variation of nonlinear coefficient term: a coefficient term vary with stress ratio in consolidation process; b variation of nonlinear terms with the stress ratio and C c /C kh(or kv) The above verifications prove that the current consolidation model based on spectral method can improve the prediction significantly especially at shallow layers where the vertical drainage can contribute considerably to the overall soil consolidation. The proposed solution is suitable for analyzing vertical and radial consolidation to capture more realistic conditions such as multilayered soils and time-dependent loading associated with nonlinear behaviors of compressibility and permeability.

Assessment of past and current nonlinear consolidation solutions
As discussed earlier, the simplified analytical solutions for nonlinear consolidation can be obtained based on certain assumptions for simplicity. While for previous models based on Assumption (1) (i.e., assuming that and Assumption (2) (i.e., assuming that r 0 r 0 are mainly determined by the ratios B h (or B v ) (i.e., -C c /C kh or -C c /C kv ) and r 0 r 0 0 . It can be seen from Fig. 10a that when the compression index (C c ) is not equal to the permeability indexes (C kv or C kh ), the nonlinear coefficient term changes significantly during the consolidation process (i.e., r 0 r 0 0 changes from 1 to r 0 0 þ q max À Á r 0 0 ). Moreover, Fig. 10b indicates that the nonlinear coefficient term changes more apparently with the increase in the effective stress ratio when C c /C kh(or kv) is less than 1. For example, the coefficient term increases sharply toward 5 when C c / C kh(or kv) = 0.5 and r 0 0 þ q max À Á r 0 0 ¼ 25: Therefore, in this section, the consolidation responses based on the Simplified Methods A and B have been obtained using the average value of . r 0 =r 0 appeared in Eq. (37), and compared with those using the proposed solution. Since the main influencing factors of the nonlinear coefficient terms are C c /C kh(or kv) and q max r 0 0 , the effects of these two ratios are investigated through the parametric study. The well resistance is neglected, and the imposed drainage condition is the PTIB (impervious bottom) with an instantaneous loading. The single layer normally consolidated soil (r 0 0 ¼ r 0 p ) is considered isotropic (C v0 = C h0 , C k = C kh = C kv ). The soil properties based on the Moruya clay (New South Wales) were assumed as follows: (i) the soil properties: C v0 = C h0 = 1.2 9 10 -3 m 2 /day, r 0 0 = 20 kPa, C c = 0.3, C k = C kh = C kv = 0.45, e 0 = 1; (ii) the permeability ratio (k h /k s ) = 1.5; and (iii) the geometrical parameters of drains: r e = 0.5 m, r s-= 0.222 m, r w = 0.074 m, n = r e /r w = 6.79, s = r s /r w-= 3.02, H = 5 m, l = 1.718. Note that series terms in relation to N = 50 were used in the analysis.

Effect of the ratio between the compression and permeability indices (C c /C k )
To study the impact of the ratio of C c /C k (x), in the range of 0.5-2 was adopted in the analysis according to Berry and Wilkinson [5], the load increment ratio q max . r 0 p was set as 5. Figure 11 (where T h = C h0 t/d e 2 ) shows the comparison between the proposed and simplified solutions for different x. Apparently, x has a great impact on the consolidation rate. It shows that given the same soil parameters and load conditions, the greater the value of x, the smaller the consolidation rate. This is because the consolidation coefficient decreases as x increases, as shown in Eqs. (7) and (8). It can also be seen that when x is greater than 1 (black and red lines), the results of the two simplified solutions are quite different from the results by the proposed method. This is because the consolidation coefficients of Simplified Methods are greater than the varying consolidation coefficient adopted in the proposed solution in the early stage, and smaller in the later stage, as shown in Fig. 11c (take x = 1.5 as an example). This results in the consolidation rate being lower in the early stage and larger in the later stage for both Simplified Methods. When x is less than 1 (green lines), an opposite trend is observed. When x is equal to 1 (blue line), i.e., the nonlinear terms (i.e., caused by the magnitude of the difference between the average consolidation coefficients adopted by Simplified Methods and the variable consolidation coefficient used in the proposed method, which can be seen from Fig. 11c and d.

Effect of load increment ratio
The load increment ratio q max =r 0 0 (R) is related to the applied preloading and the in-situ initial stress. The greater the load (q max ) or the smaller of in-situ initial stress (r 0 0 ), the greater the load increment ratio R. To study the impact of load increment ratio (R) in the range of 1-10, four different values of R were selected under the two cases of x = 1.5 and x = 0.5: (i) R = 1; (ii) R = 4; (iii) R = 7; and (iv) R = 10 ( Figs. 12 and 13).
For x = 1.5, the increase in load increment ratio reduces the consolidation rate based on EPWP (i.e., U p ), as shown in Fig. 12. In contrast, the consolidation rate increases as R increases when x = 0.5 (see Fig. 13). This is because the consolidation coefficient decreases as R increases when x [ 1, and increases with the increase of R when x \ 1, as shown in Fig. 10. It can also be seen that the results from the two simplified solutions are quite different from those by the proposed solution. When R is small (i.e., R = 1), the differences in the computational results between the simplified and the proposed solutions are relatively small (the largest difference given by both simplified solutions is less than 3.5%), as shown in Figs. 12 and 13. As the load increment ratio R increases, the deviation of the simplified solutions gradually becomes significant. When x = 1.5, Simplified Method A has a relatively small deviation in the later stage of consolidation (T h [ 0.5), and the Simplified Method B has a relatively small deviation in the early stage of consolidation (T h \ 0.1). When x = 0.5, the results of Simplified Methods A and B are relatively close, but overall, Simplified Method B has a smaller deviation. For x = 1.5 and x = 0.5, the biggest difference between the simplified solutions and proposed solution reaches 12.5% and 11.0%, respectively, when R = 10.

Applicability of the simplified solutions
The above analysis indicates that the simplified solutions can cause noticeable deviations in the predicted results depending on the magnitude of x and R. The simplified solutions must be applied in an appropriate range to maintain their prediction's accuracy. For this purpose, the typical values of C c /C k for soil in the range of 0.5-2 were used and the range of load increment ratio q max r 0 0 was selected within 0.1-10. Figure 14 shows the maximum deviations in the degrees of consolidation between the proposed solution and the Simplified Methods A and B for different values of x and R. Obviously, the deviation of both Simplified Methods A and B increases with the increase of R and |x -1|. In this study, the deviations originated by Simplified Methods A and B both reach the maximum values, i.e., 20.1% and 28.9%, respectively, when x = 2 and R = 10. If 5% error is taken as the acceptable threshold considering the deviation in predicted results, when considering the degree of consolidation based on settlement (i.e., U s ), Simplified Methods A and B can only satisfy this requirement if the following conditions are met: Combining the above conditions for a general case, it can be concluded that both Simplified Methods A and B can provide acceptable predictions below 5% error when either the load increment ratio is relatively low (R \ 2) or the compression index is close to the permeability index (0.75 \ x \ 1.10). It is noteworthily that the assumption for the smear zone would affect the value of the dimensionless parameter l. However, since the difference between the simplified solutions and the proposed solution is essentially the determination of nonlinear term , the value of l has a slight influence on the deviation from accuracy when adopting simplified solutions by additional computational verification. In this regard, the assumption for the smear zone would not change the related conclusions to any significant extent.

Model limitations
Although the proposed model can predict the nonlinear consolidation of stratified soil induced by vertical drains, it still has some limitations due to some assumptions made for facilitating the mathematical formulations and solutions. Some of these limitations are listed below: (a) The spectral-Galerkin method solution can lead to oscillations when the problem is represented by a discontinuous function; these oscillations are known as Gibbs phenomenon [67]. Therefore, more series terms are required when modeling sharp changes in the pore pressure profile. (b) The constitutive relationship associated with preloading removal has not been considered in this study.

Conclusions
In this paper, a novel approach was proposed where the spectral method was used to analyze the nonlinear consolidation of multilayered soil with coupled vertical-radial drainage. The logarithmic compressibility and permeability model (e-log r' and e-log k) was adopted to describe the nonlinear relationships. Conclusions can be drawn as follows: (1) The proposed method can capture well the nonlinear characteristics in consolidation behavior of different soil layers with time and depth. The application of this method to existing laboratory and field data in comparison with other analytical solutions verified the feasibility and accuracy of the proposed model. For the case study, the difference between the predicted and measured settlement at 400 days significantly decreased from 5.02% (i.e., 87 mm) by the previous models to 0.58% (10 mm) by the proposed model.
(2) The value of x (C c /C k ) had a great impact on the consolidation rate, i.e., the greater the value of x, the smaller the consolidation rate. Increasing the load increment ratio (R = q max r 0 0 ) and the deviation of the ratio x from unity (i.e., |x -1|) can lead to a larger deviation of both Simplified Methods A and B. The rate of strain can be expressed as: It is assumed that the flow rate in the unit cell is equal to the rate of change in the volume of the soil mass, then the continuity equation can be expressed by: where k is the radial permeability coefficient, k = k s and k h inside and outside the smear zone, respectively. The average EPWP in the soil cylinder at depth Z is calculated from the following algebraic expression: By substituting Eq. (27) into (28), the following equation expressed by the average EPWP can be obtained: where l is the dimensionless parameter, which is computed based on the variation of soil permeability within the smear zone and the radial geometry of the drain. Based on Eqs.
(2)-(5), (27) and assumptions (a)-(e), the governing equation can be expressed as: A.2: Solutions by using the spectral method By substituting Eq. (12) into Eq. (9) and using the spectral-Galerkin method, the governing differential equations can be rewritten as: By using the method of variation of parameters, the solution to the non-homogeneous Eq. (30) can be found by: To present the explicit matrix element expressions for C, W and I in a concise manner, some shorthand notations are adopted as shown below: Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.

Data availability
The data used to support the findings of this study are available from the corresponding author upon request.