Relevant research (Chen et al., 2020a; Chen et al., 2020b) has showed that the drainage volume in the middle and late stages of consolidation that engineers were concerned about was less. In order to achieve the purpose of drainage and resource conservation, scholars (Fan and Mei, 2016) proposed the method of distributing drainage bodies. The corresponding schematic diagram is shown in Figure 1. Among them, Figure 1A is a model diagram of distributed drainage boundary, and Figure 1B is a characteristic element obtained based on symmetry.

In Figure 1A, the drainage material is arranged on the top surface of the foundation at equal intervals. Among them, the width of the paved area is 2L, and the distance between adjacent paved areas (that is, the unpaved area) is 2D. In this way of laying the drainage body, the pore water in the foundation first flows horizontally from the undistributed area to the distributing area, and then flows vertically within the distributing area, and finally drains out of the foundation. In essence, the problem of ground consolidation under the distributed drainage boundary (Fan and Mei, 2016) can be attributed to the plane strain problem.

When establishing the MPSC model, the following assumptions are introduced here:

a.The soil is always saturated, and soil particles and pore water will not be compressed.

When studying the consolidation process of distributed drainage boundary foundations, scholars (Yao et al., 2019; Chen et al., 2020b) usually choose the characteristic element (i.e., Figure 1B) for modeling. Among them, the width of the characteristic element includes two parts: a paved area with a width of L and a unpaved area with a width of D. Similar to the existing research, the MPSC model also uses the above method, and the schematic diagram is shown in Figure 2. As shown in Figure 2A, the thickness of the soil layer is H ₀, and it is always in saturated state (i.e., the height of the static water surface is always H _w). At the initial moment (see Figure 2A), the soil layer is divided into R i × R j elements. Among them, the horizontal direction is R i elements of equal width (d); the vertical direction is R j elements of equal thickness (L ₀). In the analysis, the center point of the element is used as the research object. In order to facilitate the calculation of the position of the corresponding element, z coordinate system with the positive direction upward is established here, and the coordinates of the center point of the element (i.e., z i , j t ) and the coordinates of the right vertex of the element (i.e., z c , i , j t ) are given. Among them, z c , i , j t is related to the cross-sectional area of the element (i.e., A i , j t ), and the corresponding relationship is shown in Eq. 1. Moreover, z i , j t and the thickness of the consolidated element (i.e., L i , j t ) are related to z c , i , j t , and the corresponding relationship is shown in Eqs 2, 3. At the initial moment (see Figure 2A), the soil layer has been consolidated and stabilized under the original load q ₀ and its own weight. When the new load Δ q is applied, the soil layer enters a new consolidation stage (see Figure 2B). During the consolidation process, the deformation only occurs in the vertical direction, and the width of element in the horizontal direction is always the same. Regarding the drainage surface of the foundation, it is assumed that only the top surface is drained, and it can be described by a distributed drainage boundary. That is to say, the paved area is drainage area, and the corresponding hydraulic gradient is not equal to 0; the unpaved area is undrained area, and the corresponding hydraulic gradient is 0. z c , i , j t = z c , i − 1 , j t + z c , i , j − 1 t − z c , i − 1 , j − 1 t + A i , j t − z c , i − 1 , j t − z c , i − 1 , j − 1 t d d / 2 (1) z i , j t = z c , i − 1 , j t + z c , i − 1 , j − 1 t + z c , i , j t + z c , i , j − 1 t 4 (2) L i , j t = z c , i − 1 , j t − z c , i − 1 , j − 1 t + z c , i , j t − z c , i , j − 1 t 2 (3)

Regarding the plane-strain consolidation problem, (Deng and Zhou, 2016a; Deng and Zhou, 2016b) carried out corresponding research based on the piecewise-linear method. With reference to the research of Deng et al. (2016a), Deng et al. (2016b), the total stress of the corresponding element can be expressed as: σ i , j t = H w − z c , i − 1 , R j t + z c , i , R j t 2 γ w + q 0 + Δ q + A i , j t γ i , j t 2 d + ∑ n = j + 1 R j A i , n t γ i , n t d (4) Where γ w represents the specific gravity of water; γ i , j t = G s + e i , j t ) γ w / 1 + e i , j t represents the saturated weight for element (i, j). The parameter G s is the specific gravity of soil particles. Based on previous research (Fox et al., 1997) on piecewise-linear consolidation methods, the value of G s can be defined as follows: when the self-weight of the soil is not considered, G s = 1; when the self-weight of the soil is considered, the magnitude of G s depends on the category of the soil.

In MPSC model, it is assumed that the deformation of the soil has rheological characteristics, and the UH constitutive model considering the time effect (Yao et al., 2013; Hu and Yao, 2015) can be used to describe the above-mentioned characteristics of the soil. Regarding the UH constitutive model considering the time effect (Yao et al., 2013; Hu and Yao, 2015), existing studies (Liu et al., 2020) have introduced it into the one-dimensional consolidation process for corresponding discussion. Among them, the relationship between effective stress and void ratio can be expressed by Eqs 5, 6). Regarding the seepage process, it is considered that the seepage process is isotropic. And, the process conforms to Darcy’s law. σ ′ i , j t + Δ t = σ ′ i , j t 1 − e i , j t − e i , j t + Δ t ln ⁡ 10 − C α R i , j t α Δ t C s + C c − C s M 4 M f 4 , σ ′ i , j t + Δ t > σ ′ i , j t (5) σ ′ i , j t + Δ t = C s σ ′ i , j t C s − e i , j t − e i , j t + Δ t ln ⁡ 10 − C α R i , j t α Δ t , σ ′ i , j t + Δ t ≤ σ ′ i , j t (6) Where C c , C s and C α are the compression index, the swelling index and the secondary consolidation coefficient, respectively; Δ t is the time increment; M、 M f are the stress ratio at the critical state and the potential failure stress ratio, respectively. R i , j t is the over-consolidation parameter of the corresponding element at time t. The corresponding expressions for the above parameters are as follows: M = 6 ⁡ sin ⁡ φ / 3 − sin ⁡ φ (7) M f = 6 χ / R 1 + χ / R − χ / R (8) R = σ ′ / σ 0 ′ ⁡ exp − ln ⁡ 10 ε v p 1 + e 0 / C c − C s (9) Where χ = M 2 / 12 3 − M ; φ is the effective internal friction angle of the soil; σ 0 ′ is the initial effective stress adapted to the initial void ratio e ₀; It is worth noting that R ₀ is the initial value of the over-consolidation parameter R, and it has a reciprocal relationship with the over-consolidation ratio OCR; ε v p is plastic deformation.

After solving the total stress and effective stress, the corresponding pore pressure can be expressed as: u i , j t = σ i , j t − σ ′ i , j t (10)

For the seepage process, Figure 3 shows a schematic diagram of the seepage flow. It can be seen from Figure 3 that the seepage process is divided into vertical seepage and horizontal seepage. First, the vertical seepage process is introduced here. Among them, the vertical equivalent permeability coefficient can be expressed as: K z s , i , j t = k z , i , j + 1 t k z , i , j t L i , j + 1 t + L i , j t L i , j + 1 t k z , i , j t + L i , j t k z , i , j + 1 t (11) Where k z , i , j t is the vertical permeability coefficient, and its relation to void ratio can be described by e-lgk (Liu et al., 2020). e = e 0 − C k lg k 0 / k (12) Where C k is the permeability index, which can reflect the speed of change of the permeability coefficient with void ratio.

For the upper and lower boundaries, the equivalent permeability coefficients are defined as k z s , i , R j t = k v , i , R j t and k z s , i , 0 t = k z , i , 1 t .

In addition to the permeability coefficient, the vertical hydraulic gradient is i z , i , j t = h i , j + 1 t − h i , j t z i , j + 1 t − z i , j t (13) Where h i , j t is the head of element (i, j) at time t, which is composed of two parts, namely the head generated u i , j t / γ w by the pore pressure and the position head z i , j t .

In the paper, the boundary of the top surface of the soil layer can be described by distributed drainage boundary conditions, which is composed of paved area and unpaved area. In unpaved area, the top boundary of the soil layer is impervious, that is, the vertical component of the hydraulic gradient is 0. In paved area, the vertical component of the corresponding hydraulic gradient can be expressed by Eq. 14, namely: i z , i , R j t = 2 H w − h i , R j t L i , R j t ⁡ sin ⁡ θ i , R j t (14)

For horizontal seepage, the corresponding hydraulic gradient and equivalent permeability coefficient can be expressed as: i x , i , j t = h i + 1 , j t − h i , j t d 2 + z i + 1 , j t − z i , j t 2 (15) K x s , i , j t = 2 k x , i , j t k x , i + 1 , j t k x , i + 1 , j t + k x , i , j t (16) Where k x , i , j t is horizontal permeability coefficient, which is equal to the vertical permeability coefficient k z , i , j t .

Then, the flow at the boundary of the element in the corresponding direction can be expressed as: q x , i , j t = − K x s , i , j t i x , i , j t z c , i , j t − z c , i , j − 1 t sin ⁡ θ i , j t (17) q z , i , j t = − K z s , i , j t i z , i , j t d (18)

After solving the corresponding flow, related indicators (i.e., the area A i , j t + Δ t , void ratio e i , j t + Δ t , top surface settlement S i t + Δ t , Average settlement S a v g t + Δ t ) of the corresponding element at time t + Δ t can be expressed as: A i , j t + Δ t = A i , j t + q x , i − 1 , j t + q z , i , j − 1 t − q x , i , j t − q z , i , j t Δ t (19) e i , j t + Δ t = A i , j t + Δ t 1 + e 0 , i , j A 0 , i − 1 (20) S i t + Δ t = H 0 − z c , i , R j t + Δ t (21) S a v g t + Δ t = ∑ i = 1 R i S i t + Δ t R i (22) Where A 0 , j is the cell cross-sectional area at the initial moment; e 0 , i , j is the void ratio of the corresponding element at the initial moment, which is judged according to the initial stress.

In order to observe the dissipation of pore pressure in the foundation, the calculation formula of excess pore pressure is given here, namely: u e x , i , j t = u i , j t − u 0 , i , j (23) Where u 0 , i , j is the pore pressure of the corresponding element at the initial moment.

Then, the average excess pore pressure in the foundation can be expressed as: u e x , a v g t = ∑ i = 1 R i ∑ j = 1 R j d L i , j t u e x , i , j t L + D H 0 (24)

In order to simulate the rheology of soil, many scholars have proposed corresponding constitutive relations. Among them, Yao et al. (2013) proposed the UH model considering the time effect. Later, Hu and Yao, (2015) introduced its one-dimensional form into the classic Terzaghi consolidation theory and studied the effect of rheology on the consolidation process. In order to verify the validity of the MPSC model in this paper, the case provided by Hu and Yao, (2015) is used for comparison. When making comparison, the MPSC model needs to be degenerated into a one-dimensional form. The specific method is: the number of horizontal units is 1, and the top drainage boundary is a completely permeable boundary. What’s more, the values of other parameters are: the soil thickness H ₀=1m, H _w=1m, G s =1.0, q ₀=10kPa, e ₀=0.53, Δ q =90 kPa. Regarding the constitutive model, the corresponding parameter value is as follows: C c =0.0217, C s =0.0131, M =1.112, R ₀=0.95. For the secondary consolidation coefficient, the corresponding values are 0.5 C α , C α and 2 C α . Among them, C α =0.0108. Regarding the seepage process, Darcy’s law is considered to be applicable, and the permeability coefficient k=3.63 × 10^-7 m/min and remains unchanged.

Figure 4 shows the comparison results of the average consolidation degree U _P obtained by Hu and Yao, (2015) and MPSC model in this paper. By comparing the data in Figure 4, it is found that the maximum error between the calculated results of the two methods is 2.6%. That is to say, the calculation results in this paper can better match the results of Hu and Yao, (2015). This also indicates that the calculation model in this paper is effective for the rheological consolidation calculation process.

In order to optimize the layout of horizontal drainage bodies, scholars (Fan and Mei, 2016) proposed the distributed drainage boundary. Existing studies (Yao et al., 2019; Chen et al. (2020a), Chen et al. (2020b) have shown that the boundary can achieve the purpose of drainage and resource utilization. In previous studies, Chen et al. (2020b) used analytical solutions to explore the influence of the width-thickness ratio Q on the pore pressure dissipation process under the linear elastic constitutive relationship. Among them, time is expressed in a dimensionless way, and the corresponding expression is T v = C v t H 0 2 = k 1 + e 0 a γ w H 0 2 t (25) Where a is the compression factor; e ₀ is the initial void ratio.

In order to verify the correctness of the MPSC model in this paper, the above case is selected for comparison. When making comparison, the flow process of water in soil pores conforms to Darcy’s law, and the UH model considering the time effect in MPSC model needs to replaced by the elastic constitutive relationship (i.e., Δ e = a Δ p ) as Chen et al. (2020b). Additionally, the other conditions set in this paper in the comparison are consistent with those used by Chen et al. (2020b). The comparison results are shown in Figure 5, and the maximum error between the calculated results of the two methods is 3.4%. This indicates that the MPSC model is capable of simulating soil consolidation problems under distributed drainage boundary conditions.

In the study of distributed drainage boundary, the pave rate W is a key parameter. Generally, the pave rate refers to the ratio of the width of the paving drainage body to the total width of the soil layer. In this case, the pave rate W is taken as 0.1, 0.5 and 0.8 in sequence. Figure 6 shows the effect of the pave rate W on the average excess pore pressure u e x , a v g . It can be found from this figure that at the initial stage of consolidation (such T _v=0.001), the average excess pore pressure u e x , a v g is greater than the applied load (i.e., new load Δ q ). This is a phenomenon that has not been discovered by existing studies (Fan and Mei, 2016; Yao et al., 2019; Chen et al., 2020b) under distributed drainage boundary. Regarding the cause of this phenomenon, most one-dimensional studies (Hu and Yao, 2015; Liu et al., 2020) attribute it to rheological effects. Besides, when the boundary is a distributed drainage boundary (that is, the pave rate W is less than 1.0), the pore pressure curve is always above the full-distribution drainage boundary. This shows that the pave rate can delay the dissipation of pore pressure in the foundation. Moreover, the degree of delay will increase as the pave rate decreases. For example, at T _v = 0.4, the average excess pore pressure u e x , a v g with the pave rates of 0.1, 0.5 and 0.8 are 95.01 kPa, 68.18 kPa, and 46.28 kPa, respectively. However, when the pave rate is 0.8, the pore pressure curve will gradually approach the pore pressure curve at the boundary of the completely arranged drainage body. This also shows that reasonable laying of horizontal drainage materials can achieve the effect of both drainage and resource conservation.

Figure 7 shows the influence of the pave rate W on the average settlement S a v g . On the whole, the smaller the pave rate, the smaller the settlement at the same time. For instance, at T _v = 0.2, when the pave rate is 0.1, 0.5 and 0.8, the values of S a v g are 0.072 m, 0.156 m, 0.226 m, respectively. Moreover, they are respectively 0.28 times, 0.61 times, and 0.88 times of the corresponding settlement under the full-distributed drainage boundary. In other words, to reach the same value of settlement, the greater the pave rate, the shorter the time consumed. Taking S a v g = 0.3 m as an example, the time required for the pave rate equal to 0.8 when the settlement value is reached is only 0.26 times the pave rate equal to 0.1.

In addition to the pave rate, the thickness-width ratio Q is also an important parameter in the study of distributed drainage boundaries. In this case, the thickness-to-width ratio is 1, 2, and 6 respectively to carry out the corresponding research. Figure 8 shows the effect of Q on u e x , a v g . Similar to the above results, after considering the effect of thickness-to-width ratio, the average pore pressure dissipation curve is also located above the pore pressure curve of the full-distributed drainage boundary. This also means that the thickness-to-width ratio can also delay the dissipation of pore pressure. Moreover, the smaller the thickness-to-width ratio, the more obvious the above-mentioned delay effect. For example, when T _v= 0.5, the average excess pore pressure when the thickness-to-width ratio is 1 is still 1.42 times the corresponding pore pressure when the value of Q is 6. It can be seen that the influence of the thickness-to-width ratio on the dissipation of pore water pressure cannot be ignored.

Figure 9 shows the thickness-to-width ratio on the average settlement S a v g . It can be found from this figure that the ratio mainly affects the middle and late stages of consolidation. And, as the thickness-to-width ratio decreases, the value of settlement at the same time is smaller. For example, when T _v = 0.3, the settlement when the thickness-to-width ratio is equal to 1 is only 0.73 times the corresponding settlement when the thickness-to-width ratio is equal to 6. Furthermore, it can be found from the figure that the corresponding curves considering the thickness-to-width ratio are all located above the settlement curve of the full-distribution drainage boundary. In other words, the settlement process of the foundation after considering the thickness-to-width ratio will be slower than the settlement process under the full-distributed drainage boundary.

The rheological properties of cohesive soil is an important feature that distinguishes it from other soils. The existing researches (Liu et al., 2020; Zhou et al., 2020) indicate that the coefficient of secondary consolidation is an important indicator for describing the rheological characteristics of soil. Based on this, with the secondary consolidation coefficient C α as the representative, the influence of the rheology on the consolidation process under the distributed drainage boundary is analyzed here. In the analysis, the value of C α is taken as 0.006, 0.012, and 0.018 in sequence.

Figure 10 shows the effect of the secondary consolidation coefficient on the average excess pore pressure. It can be found from this figure that the greater the secondary consolidation coefficient, the greater the peak average of the excess pore pressure. For example, the peak of the pore pressure when C α = 0.018 is 1.23 times the corresponding peak of C α =0.006. The above phenomenon has been discovered by scholars (Hu and Yao, 2015; Liu et al., 2020) in one-dimensional consolidation analysis. Moreover, scholars (Hu and Yao, 2015; Liu et al., 2020) have also provided corresponding explanations, that is, the rheological characteristics of soil cause obstruction in the discharge of pore water in the soil.

Figure 11 shows the effect of secondary consolidation coefficient on the average settlement S a v g . It can be seen from this figure that the larger the secondary consolidation coefficient (that is, the stronger the rheological effect), the larger the value of settlement at the same time. For example, at T _v =1.0, with the increase of the secondary consolidation coefficient, the value of the average settlement is 0.334 m, 0.391 m, 0.443 m, respectively. From this result, the settlement when the secondary consolidation coefficient is equal to 0.018 is 1.33 times the settlement when the secondary consolidation coefficient is 0.006. It can be seen that the rheological effect has a significant influence on the settlement of the distributed drainage boundary foundation. Besides, it can also be seen from Figure 11 that the greater the secondary consolidation coefficient, the greater the slope of the end of the deformation curve, and the longer it takes for the settlement of the foundation to reach a stable state. From the above analysis, it can be seen that it is very necessary to accurately measure the secondary consolidation coefficient of soil when predicting the settlement of soil in practical engineering.

The over-consolidation of soil is an important branch of consolidation research. In the UH model considering the time effect (Yao et al., 2013; Hu and Yao, 2015; Liu et al., 2020), the initial over-consolidation parameter R ₀ reflect the initial over-consolidation state of the soil. In the study, the initial over-consolidation parameter R ₀ is 0.3, 0.6 and 0.9 in turn. Figure 12 shows the effect of the parameter R ₀ on the average excess pore pressure. From this figure, the initial over-consolidation state has significant effect on the dissipation process of the excess pore pressure in the distributed drainage boundary foundation. On the whole, the larger the value of R ₀ (that is, the weaker the initial over-consolidation state), the slower the average pore pressure will dissipate. For example, at T _v = 0.2, with the increase of R ₀, the corresponding average pore pressures are 28.62 kPa, 58.60 kPa and 88.25 kPa. In other words, when it dissipates to the same pore pressure value, the stronger the initial over-consolidation of the soil, the shorter the time it takes. For example, when the average excess pore pressure dissipates to 52 kPa, the dimensionless time taken when the initial over-consolidation parameters are 0.3, 0.6, and 0.9 are 0.10, 0.26, and 0.60, respectively.

Figure 13 shows the influence of the parameter R ₀ on average settlement. Figure 13 shows that the initial over-consolidation state of the soil has a significant impact on the settlement process of the foundation. For the same moment, the smaller the initial over-consolidation parameter R ₀, the smaller the corresponding average settlement. For example, when the dimensionless time T _v = 1, as the initial over-consolidation parameter increases, the corresponding settlement values are 0.129 m, 0.281 m, and 0.389 m. It can be seen that the settlement at this moment when the initial over-consolidation parameter is 0.9 is 3.02 times the corresponding settlement when the initial over-consolidation parameter is 0.3. In addition, it can also be seen from the figure that the stronger the degree of initial over-consolidation, the shorter the time it takes for the foundation to enter a stable state. For example, when the dimensionless time T _v is equal to 1, the foundation has reached a stable state when the initial over-consolidation parameter is equal to 0.3, but the deformation of the foundation when the initial over-consolidation parameter is equal to 0.6 and 0.9 is still in progress. This also indicates that preloading is an effective method to reduce post construction settlement.