Viscoelastic Parameter Prediction of Multi-Layered Coarse-Grained Soil with Consideration of Interface-Layer Effect

Study on viscoelastic properties of the multi-layered coarse-grained soil (CGS) is very important for safety assessment and disaster prevention of subgrade engineering. Current research work is mainly focused on the one-layered CGS and the actual pebble inclusion of irregular polyhedron is usually simplified as an ideal shape of sphere or ellipsoid. Very few studies are available for predicting viscoelastic parameters of the multi-layered CGS. In this paper, a new method is proposed to predict viscoelastic parameters of multi-layered CGS based on the homogenization method and elastic–viscoelastic corresponding principle, in which the interface-layer viscoelasticity and the actual shape of pebble inclusion are firstly taken into account. Research results show the creep deformation is decreased with the increase of the shape factor (ρ) of pebble inclusion, and the interface-layer height (h) and numbers (N). ρ is in the range of 1–1.8 and the suitable interface-layer height is 20–30% as much as the height of one-layered CGS. The tested creep curves of the multi-layered CGS agree well with the predicted ones and can prove the existence of the interface-layer (considering at least one interface-layer) and verify the validity of this new interface-layer method.


Introduction
Coarse-grained soil (CGS) consists of a pebble (with high elastic modulus) and sand-clay matrix (with rheological properties) surrounding the pebble and pore [1][2][3]. It is easily available in nature and thus widely applied in the subgrade, bridge foundation, and structure foundation engineering [4]. For improving the bearing capacity of the foundation, the CGS of different gradations needs to be artificially filled [5][6][7][8]. The slipping and embedding of different gradations CGS particles would lead to a certain height of interface-layer, i.e., the formation of multi-layered CGS [9][10][11]. Under long-term loading, the multi-layered CGS may cause uneven settlement due to its rheology properties (depending on its microscopic components and structure) and there exists a potential safety hazard for the foundation engineering [12]. Therefore, it is very important to determine the macroscopic rheological parameters of multi-layered CGS with consideration of the interface-layer for predicting its long-term deformation.
Currently, the macroscopic rheological parameters of materials are mainly studied by two theoretical methods: periodic cell method and homogenization method. Since the periodic cell method needs a representative volume unit of a periodic array, it is only used for most metallic alloy materials [13], ceramic material [14], and a few nonmetallic materials with periodic microstructures (e.g., concrete [15]). The homogenization method is widely applied for the geotechnical materials Figure 1 shows a diagrammatic sketch of three-layered CGS composite structures with solid phase (pebble, sand-clay matrix) and gas phase (pore) of different gradations, where H is one-layered CGS height, and h 1 and h 2 are interface-layer height. For determining its macroscopic viscoelastic parameters, its macroscopic elastic parameters need to be firstly obtained based on homogenization principle. This includes four steps: (1) consider the solid phase (sand and clay as matrix, pebble as inclusion) for homogenizing the one-layer CGS without consideration of gas phase (pore), in order to obtain the effective elastic parameters (bulk and shear modulus) of the solid phase, (2) consider both the solid phase (as matrix) and the gas phase (pore as inclusion) for homogenization in order to obtain the elastic parameters (bulk and shear modulus) of the one-layered CGS, (3) consider the interface-layer as two one-layers with different gradations of the solid phase and the gas phase for homogenization, in order to obtain the elastic parameters (bulk and shear modulus) of the interface-layer, (4) superimpose the results of steps (2) and (3) to obtain the effective elastic modulus of multi-layered CGS. Then, apply the elastic-viscoelastic corresponding principle to obtain the viscoelastic parameters of multi-layered CGS, details as follows. Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 17

Viscoelastic Parameter Model of One-Layered CGS
Effective elastic parameters of the one-layered CGS can be calculated by those of the solid phase with consideration of the gas phase (pore). Figure 2 shows the calculation model of the solid phase for one-layered CGS, where Ωp is the region of pebble, Ωm is the region of sand-clay matrix and Ωeq is the region of equivalent media. The pebble is first simplified as isotropic sphere inclusion and completely surrounded by the isotropic sand-clay matrix. According to the generalized self-consistent model, the effective bulk modulus (Ks) and shear modulus (Gs) can be obtained [19,36].
where Kp and Km, Gp and Gm are the bulk modulus and shear modulus of pebble and sand-clay matrix, p0 is the volume fraction of pebble, A1, B1, C1 are undetermined coefficients. It is worth mentioning that the pebble is actually an irregular polyhedron rather than the regular sphere. Under the same volume condition, the sphere has the smallest surface area while the tetrahedron (the simplest convex polyhedrons) has the largest surface area. Since the contacted surface area of pebble with sand-clay matrix affects the mechanical properties of the solid phase greatly, the volume of sphere pebble must be increased in order to match the same surface area of the actual irregular convex polyhedron. Define a shape factor ρ to describe the amplification effect of actual polyhedron on the volume fraction of sphere pebble. For the same surface area, the sphere

Viscoelastic Parameter Model of One-Layered CGS
Effective elastic parameters of the one-layered CGS can be calculated by those of the solid phase with consideration of the gas phase (pore). Figure 2 shows the calculation model of the solid phase for one-layered CGS, where Ω p is the region of pebble, Ω m is the region of sand-clay matrix and Ω eq is the region of equivalent media. The pebble is first simplified as isotropic sphere inclusion and completely surrounded by the isotropic sand-clay matrix.

Viscoelastic Parameter Model of One-Layered CGS
Effective elastic parameters of the one-layered CGS can be calculated by those of the solid phase with consideration of the gas phase (pore). Figure 2 shows the calculation model of the solid phase for one-layered CGS, where Ωp is the region of pebble, Ωm is the region of sand-clay matrix and Ωeq is the region of equivalent media. The pebble is first simplified as isotropic sphere inclusion and completely surrounded by the isotropic sand-clay matrix. According to the generalized self-consistent model, the effective bulk modulus (Ks) and shear modulus (Gs) can be obtained [19,36].
where Kp and Km, Gp and Gm are the bulk modulus and shear modulus of pebble and sand-clay matrix, p0 is the volume fraction of pebble, A1, B1, C1 are undetermined coefficients.
It is worth mentioning that the pebble is actually an irregular polyhedron rather than the regular sphere. Under the same volume condition, the sphere has the smallest surface area while the tetrahedron (the simplest convex polyhedrons) has the largest surface area. Since the contacted surface area of pebble with sand-clay matrix affects the mechanical properties of the solid phase greatly, the volume of sphere pebble must be increased in order to match the same surface area of the actual irregular convex polyhedron. Define a shape factor ρ to describe the amplification effect of actual polyhedron on the volume fraction of sphere pebble. For the same surface area, the sphere According to the generalized self-consistent model, the effective bulk modulus (K s ) and shear modulus (G s ) can be obtained [19,36].
where K p and K m , G p and G m are the bulk modulus and shear modulus of pebble and sand-clay matrix, p 0 is the volume fraction of pebble, A 1 , B 1 , C 1 are undetermined coefficients. It is worth mentioning that the pebble is actually an irregular polyhedron rather than the regular sphere. Under the same volume condition, the sphere has the smallest surface area while the tetrahedron (the simplest convex polyhedrons) has the largest surface area. Since the contacted surface area of pebble with sand-clay matrix affects the mechanical properties of the solid phase greatly, the volume of sphere pebble must be increased in order to match the same surface area of the actual irregular convex polyhedron. Define a shape factor ρ to describe the amplification effect of actual polyhedron on the Appl. Sci. 2020, 10, 8879 4 of 16 volume fraction of sphere pebble. For the same surface area, the sphere radius of r and the tetrahedron length of l have the relationship as Equation (2): where ρ max is the maximum value of shape factor ρ.
For the actually irregular convex polyhedron of pebble, ρ = 1-1.8. Substituting ρ into Equation (1) leads to the effective bulk and shear modulus of the solid phase, as Equation (3): According to Equation (4): the effective elastic modulus (E s ) of the solid phase is shown as Equation (5): For calculating the effective elastic modulus of the one-layered CGS, gas phase (pore) must be considered. Figure 3 shows the calculation model of the one-layered CGS, i.e., the hollow sphere model, where Ω p1 is the region of gas phase (pore) as inclusion, Ω m1 is the region of solid phase (sand-clay-pebble matrix) and Ω eq1 is the region of equivalent media. radius of r and the tetrahedron length of l have the relationship as Equation (2) (2) where ρmax is the maximum value of shape factor ρ.
For the actually irregular convex polyhedron of pebble, ρ = 1-1.8. Substituting ρ into Equation (1) leads to the effective bulk and shear modulus of the solid phase, as Equation (3): According to Equation (4): the effective elastic modulus (Es) of the solid phase is shown as Equation (5) For calculating the effective elastic modulus of the one-layered CGS, gas phase (pore) must be considered. Figure 3 shows the calculation model of the one-layered CGS, i.e., the hollow sphere model, where Ωp1 is the region of gas phase (pore) as inclusion, Ωm1 is the region of solid phase (sandclay-pebble matrix) and Ωeq1 is the region of equivalent media. According to the hollow sphere model, the effective bulk modulus (Ko) and shear modulus (Go) of one-layer CGS with two phases of solid and gas can be obtained [16,36] as Equation (6): (1 ) / (4 3 ) where p1 is pore volume fraction, A2, B2, C2 are undetermined coefficients. Similarly, the elastic modulus (Eo) is shown as Equation (7). According to the hollow sphere model, the effective bulk modulus (K o ) and shear modulus (G o ) of one-layer CGS with two phases of solid and gas can be obtained [16,36] as Equation (6): where p 1 is pore volume fraction, A 2 , B 2 , C 2 are undetermined coefficients. Similarly, the elastic modulus (E o ) is shown as Equation (7).
According to elastic-viscoelastic corresponding principle (the boundary value problem of linear viscoelastic body has the same solution form as that of elastic body in the Laplace space. Therefore, the viscoelastic problem is firstly transformed into Laplace space to calculate its elastic solution and then Appl. Sci. 2020, 10, 8879 5 of 16 inversely transformed into Euler space to determine its viscoelastic parameters) [37,38], the viscoelastic parameter of one-layered CGS can be obtained as Equation (8).
where the superscript T means Laplace transform, s is the Laplace transform parameter.
The viscoelastic parameter of one-layered CGS is shown as Equation (9): Figure 4 shows the calculation model of the interface-layer, where Ω p2 is the region of gas phase (pore) as inclusion, Ω m2 and Ω m3 are the regions of two-layer solid phases and Ω eq2 is the region of equivalent media. Since Ω m2 and Ω m3 are known from the above calculation (Equations (3)-(5)), the interface-layer model can be equivalent to superposition of the hollow sphere model (Ω p2 + Ω m2 ) and generalized self-consistent model (Ω m3 + Ω p2 + Ω m2 ).

Viscoelastic Parameter Model of the Interface-Layer
According to elastic-viscoelastic corresponding principle (the boundary value problem of linear viscoelastic body has the same solution form as that of elastic body in the Laplace space. Therefore, the viscoelastic problem is firstly transformed into Laplace space to calculate its elastic solution and then inversely transformed into Euler space to determine its viscoelastic parameters) [37,38], the viscoelastic parameter of one-layered CGS can be obtained as Equation (8).
where the superscript T means Laplace transform, s is the Laplace transform parameter.
The viscoelastic parameter of one-layered CGS is shown as Equation (9): Figure 4 shows the calculation model of the interface-layer, where Ωp2 is the region of gas phase (pore) as inclusion, Ωm2 and Ωm3 are the regions of two-layer solid phases and Ωeq2 is the region of equivalent media. Since Ωm2 and Ωm3 are known from the above calculation (Equations (3)-(5)), the interface-layer model can be equivalent to superposition of the hollow sphere model (Ωp2 + Ωm2) and generalized self-consistent model (Ωm3 + Ωp2 + Ωm2). According to the hollow sphere model, the equivalent bulk modulus (Ke) and shear modulus (Ge) of Ωp2 − Ωm2 region can be obtained as Equation (10): (1

Viscoelastic Parameter Model of the Interface-Layer
where Ke and Ge are the one-layered CGS effective bulk modulus and shear modulus, p2 is the volume fraction of Ωp2 − Ωm2 region, A3, B3, C3 are undetermined coefficients. According to the generalized self-consistent model, the equivalent bulk modulus (Ki) and shear modulus (Gi) of interface-layer ((Ωp2 + Ωm2) − Ωm3 region) can be calculated as Equation (11): According to the hollow sphere model, the equivalent bulk modulus (K e ) and shear modulus (G e ) of Ω p2 − Ω m2 region can be obtained as Equation (10): where K e and G e are the one-layered CGS effective bulk modulus and shear modulus, p 2 is the volume fraction of Ω p2 − Ω m2 region, A 3 , B 3 , C 3 are undetermined coefficients. According to the generalized self-consistent model, the equivalent bulk modulus (K i ) and shear modulus (G i ) of interface-layer ((Ω p2 + Ω m2 ) − Ω m3 region) can be calculated as Equation (11): where p 3 is the volume fraction of pore and A 4 , B 4 , C 4 are undetermined coefficients. Similarly, effective elastic parameters of the interfacial layer (E i ) is shown as Equation (12): Similarly, the viscoelastic parameter (J i ) in Laplace space is shown as Equation (13): Appl. Sci. 2020, 10, 8879 6 of 16 The viscoelastic parameter of interface-layer is shown as Equation (14):

Viscoelastic Parameter Formulae of the Multi-Layered CGS
Since the multi-layered CGS consists of several one-layered CGS and interface-layers, its elastic modulus (E) can be obtained by the series formula as Equation (15): where E u , E i , E l are the elastic modulus of the upper layer, interface-layer, and lower layer, respectively. According to the elastic-viscoelastic corresponding principle, there is the relationship between the elastic parameters (E T ) and viscoelastic parameters (J T ) (Equation (16)) in the Laplace space: For determination of the viscoelastic parameters J in Euler space, E T is firstly calculated by the Laplace transform of E (Equation (15)) and then J T is obtained by Equation (16). Finally, J can be calculated by the inverse Laplace transform of J T , as Equation (17):

Measuring Viscoelastic Parameters of the Sand-Clay Matrix
According to the Code for Design on Subgrade of Railway in China, TB10001-2005 [39], three different mass ratios of sand to clay (Table 1) were selected and uniformly mixed for preparation of the cylindric specimens (Φ 300 mm × 600 mm) by layered compacting method ( Figure 5). Triaxial creep tests were conducted by TAJ-2000 triaxial testing system (Figure 6), where the axial stress σ 1 was 300 kPa and confining pressure σ 3 was 200 kPa (σ 3 < σ 1 ).      Figure 7 shows the tested creep curves of different sand-clay matrices. They are all composed of three stages: instantaneous deformation (t = 0) stage, attenuation deformation stage (slope is decreased) and stable creep deformation stage (slope is constant). The higher the clay content, the larger the instantaneous deformation and the longer the duration time of the attenuation deformation stage. The clay content has a smaller effect on the stable creep deformation rate, since the elastic properties of the sand-clay matrices are mainly decided by the sand content and the viscous properties by clay content, that is, the viscous properties will play a dominant role as time continues, especially after the attenuation deformation stage, so the creep curves have an almost identical slope on stable creep deformation for viscosity of clay.  Figure 7 shows the tested creep curves of different sand-clay matrices. They are all composed of three stages: instantaneous deformation (t = 0) stage, attenuation deformation stage (slope is decreased) and stable creep deformation stage (slope is constant). The higher the clay content, the larger the instantaneous deformation and the longer the duration time of the attenuation deformation stage. The clay content has a smaller effect on the stable creep deformation rate, since the elastic properties of the sand-clay matrices are mainly decided by the sand content and the viscous properties by clay content, that is, the viscous properties will play a dominant role as time continues, especially after the attenuation deformation stage, so the creep curves have an almost identical slope on stable creep deformation for viscosity of clay. Burgers rheological model is usually adopted to describe the creep curve as Equation (18): where Jm is the parameter of Burgers model, E1 and E2 are elastic parameters, η1 and η2 are the viscosity parameters, t is time. Table 2 lists the fitting viscoelastic parameters of different sand-clay matrices by the Burgers rheological model with high precision (R 2 > 0.98). They are all decreased as the clay content increases, indicating that the viscoelastic properties of the sand-clay matrices mainly depend on the clay contents.  Burgers rheological model is usually adopted to describe the creep curve as Equation (18): where J m is the parameter of Burgers model, E 1 and E 2 are elastic parameters, η 1 and η 2 are the viscosity parameters, t is time. Table 2 lists the fitting viscoelastic parameters of different sand-clay matrices by the Burgers rheological model with high precision (R 2 > 0.98). They are all decreased as the clay content increases, indicating that the viscoelastic properties of the sand-clay matrices mainly depend on the clay contents.

Predicting Viscoelastic Parameters of Multi-Layered CGS
In the multi-layered CGS, the pebble and the sand-clay matrix are regarded as elastic and viscoelastic material, respectively. They are exploited from the local area, with the elastic modulus of pebble (Ep = 50 GPa) and density is 2.54 g/cm 3 , and Poisson's ratio of pebble (ν p = 0.25) and sand-clay matrix (ν m = 0.25). Assuming the volume fraction of pore (p 1 ) is 20% [40], the volume fraction of pebble (p 0 ) and sand-clay matrix (p 3 ) can be calculated by the mass-density formula.
Multi-layered CGS consists of several one-layered CGS and interface-layers (Figure 1). The volume fraction of one-layered CGS and the basic parameters of each component are in Table 3. Table 3. Basic parameters of one-layered CGS. Based on Table 3, viscoelastic parameters (ρ = 1) of the multi-layered CGS can be calculated by the following steps:

1.
Calculate the elastic modulus (E T , superscript T means Laplace space) of the sand-clay matrix by the elastic-viscoelastic correspondence principle as Equation (19).
where s is Laplace transform parameter.
Calculate the equivalent elastic modulus (K i , E i ) of the interface-layers by Equations (10)-(12).
Calculate the equivalent viscoelastic parameters (J) of the multi-layered CGS by Equation (16). and elastic-viscoelastic correspondence principle, as listed in Table 4.

Experimental Verification
Triaxial creep tests were conducted to verify the viscoelastic parameters of multi-layered CGS (Table 4), since the predicted creep curves could be calculated by the predicted viscoelastic parameters of multi-layered CGS for comparison of test results.
Take three-layered CGS as an example and assume that each single layer has the same heights (H) and each interface-layer (height h 1 and h 2 ) has the same proportional height (0.5 h 1 or 0.5 h 2 ) in upper and lower layers (Figure 1). Therefore, the actual heights of the top and bottom layers are H = H − 0.5 h 1 and H − 0.5 h 2 and the middle layer is H = H − 0.5h 1 − 0.5h 2 .

Triaxial Creep Test
Cylindric specimen of Φ 300 mm × 600 mm in size was adopted in the triaxial creep test. Heights of each layer were 300 mm for two-layered CGS and 200 mm for three-layered CGS. For preparation of one-layered CGS specimen, the pebble was added into the sand-clay matrix (Table 3) to uniformly mixed together by layered compacting method. Table 3 and Figure 8 show their volume fraction and tested gradation curves of pebble, sand, and clay. One-layered, two-layered, and three-layered CGS specimens were tested and their layer combinations were listed in Table 4. of each layer were 300 mm for two-layered CGS and 200 mm for three-layered CGS. For preparation of one-layered CGS specimen, the pebble was added into the sand-clay matrix (Table 3) to uniformly mixed together by layered compacting method. Table 3 and Figure 8 show their volume fraction and tested gradation curves of pebble, sand, and clay. One-layered, two-layered, and three-layered CGS specimens were tested and their layer combinations were listed in Table 4. Similar to Section 3.1, TAJ-2000 triaxial testing system was used in triaxial creep test. The axial stress σ1 was 300kPa and confining pressure σ3 was 200 kPa. Figure 9 shows some specimens after step loading (last step stress is σ1 = 500 kPa) creep tests. Similar to Section 3.1, TAJ-2000 triaxial testing system was used in triaxial creep test. The axial stress σ 1 was 300kPa and confining pressure σ 3 was 200 kPa. Figure 9 shows some specimens after step loading (last step stress is σ 1 = 500 kPa) creep tests.
(c) Similar to Section 3.1, TAJ-2000 triaxial testing system was used in triaxial creep test. The axial stress σ1 was 300kPa and confining pressure σ3 was 200 kPa. Figure 9 shows some specimens after step loading (last step stress is σ1 = 500 kPa) creep tests.  Figure 10 shows tested creep curves (ε-t curve) of the one-layered CGS as well as their predicted ones with different values of the shape factor ρ (ρ = 1~1.8, Equation (2)) for comparison. It is seen that  Figure 10 shows tested creep curves (ε-t curve) of the one-layered CGS as well as their predicted ones with different values of the shape factor ρ (ρ = 1~1.8, Equation (2)) for comparison. It is seen that all of the tested creep curves consist of three stages: instantaneous deformation, attenuation deformation, and stable creep deformation. Obviously, the higher the content (volume fraction) of pebble (L1 > L2 > L3 in Table 4), the smaller the creep deformation of the one-layered CGS (Figure 10a < Figure 10b < Figure 10c).

Creep Curve of One-Layered CGS
Appl. Sci. 2020, 10, x FOR PEER REVIEW 11 of 17 all of the tested creep curves consist of three stages: instantaneous deformation, attenuation deformation, and stable creep deformation. Obviously, the higher the content (volume fraction) of pebble (L1 > L2 > L3 in Table 4), the smaller the creep deformation of the one-layered CGS (Figure 10a < Figure 10b < Figure 10c). It is also found that all of the predicted creep curves of the one-layered CGS decreased with the increase of ρ, since the larger the value of ρ (i.e., the more irregular the shape of pebble), the larger the contacted surface area of pebble with sand-clay matrix and hence the smaller the creep deformation. Furthermore, the predicted curves of ρ = 1 (ideal sphere pebble) are all larger than the tested curves. That is because the actual irregular polyhedron of pebble with a large contacted surface area would lead to small creep deformation.
Comparison of the tested and predicted creep curves of different sand-clay matrix (L1, L2, L3), indicates they are all in good agreement at ρ = 1.35, 1.3, 1.25, respectively. This means the one-layered CGS with a higher content of pebble (L1 > L2 > L3) has a greater possibility of irregular pebble and thus a larger value of ρ.  Table 4, i.e., J-t curves. Similar to J-t curve of the one-layered CGS, they all consist of three stages: instantaneous It is also found that all of the predicted creep curves of the one-layered CGS decreased with the increase of ρ, since the larger the value of ρ (i.e., the more irregular the shape of pebble), the larger the contacted surface area of pebble with sand-clay matrix and hence the smaller the creep deformation. Furthermore, the predicted curves of ρ = 1 (ideal sphere pebble) are all larger than the tested curves. That is because the actual irregular polyhedron of pebble with a large contacted surface area would lead to small creep deformation.
Comparison of the tested and predicted creep curves of different sand-clay matrix (L1, L2, L3), indicates they are all in good agreement at ρ = 1.35, 1.3, 1.25, respectively. This means the one-layered CGS with a higher content of pebble (L1 > L2 > L3) has a greater possibility of irregular pebble and thus a larger value of ρ. Figure 11 shows viscoelastic parameter-time curves of three one-layered CGS (L1, L2, L3) and their combined interface-layers (L1/L2, L2/L3, L1/L3) which can be drawn by data in Table 4, i.e., J-t curves. Similar to J-t curve of the one-layered CGS, they all consist of three stages: instantaneous deformation, attenuation deformation, and stable creep deformation. The J of the interface-layers (e.g., L1/L2) are between those of the two combined one-layered CGS (e.g., L1 and L2) and smaller than the average values of two one-layered CGS (e.g., (L1 + L2)/2).  Figure 12 shows the tested creep strains (ε-t) of the two-layered CGS (L1 and L2, L2 and L3, L1 and L3) as well as the predicted ε-t curves with different interface-layer heights, (h = 0~0.5 H, H is the designed height of each layered CGS, H = 300 mm for the two-layered CGS, the total specimen height is 600 mm). They all have the same stages as those of the one-layered CGS: instantaneous deformation, attenuation deformation, and stable creep deformation. In addition, they all decrease  Figure 12 shows the tested creep strains (ε-t) of the two-layered CGS (L1 and L2, L2 and L3, L1 and L3) as well as the predicted ε-t curves with different interface-layer heights, (h = 0~0.5 H, H is the designed height of each layered CGS, H = 300 mm for the two-layered CGS, the total specimen height is 600 mm). They all have the same stages as those of the one-layered CGS: instantaneous deformation, attenuation deformation, and stable creep deformation. In addition, they all decrease as the interface-layer height (h) increases. That is because the interface-layer has a smaller viscoelastic parameter than the average values of the two one-layered CGS, e.g., J L1/L2 < (J L1 + J L2 )/2, the total creep strain (ε) is decreased with increase of h according to the following formula (Equation (21)), when H is given. The higher the interface-layer, the weaker the rheological behavior of the two-layered CGS.

Creep Curve of Two-Layered CGS
Appl. Sci. 2020, 10, x FOR PEER REVIEW 13 of 17 It is also found that most of the tested creep strain are between the predicted ε-t curves of 20%H and 30% H, which can prove the existence of the interface-layer and validity of the interface-layer calculation model. The suitable interface-layer height is h = (20-30%)H.  Figure 13 shows the tested creep strain of the three-layered CGS (L1, L2 and L3) as well as its predicted ε-t curves with consideration of different interface-layer numbers: N = 1 (L1 + L1/L2 + L2 + L3 or L1 + L2 + L2/L3 + L3) and N = 2 (L1 + L1/L2 + L2 + L2/L3 + L3), where each interface-layer height is h = 0.3 H (H = 200 mm for the three-layered CGS, the total specimen height is 600 mm). Similarly, they all have the same three stages as those both the one-layered and two-layered CGS: instantaneous deformation, attenuation deformation, and stable creep deformation.

Creep Curve of Three-Layered CGS
The three-layered CGS has two interface-layers. The more the interface-layers are taken into account in calculation equation (Equation (19)), the predicted creep deformation becomes smaller, since it is known from Section 4.2 that the interface-layer with some height would decrease the It is also found that most of the tested creep strain are between the predicted ε-t curves of 20%H and 30% H, which can prove the existence of the interface-layer and validity of the interface-layer calculation model. The suitable interface-layer height is h = (20-30%)H. Figure 13 shows the tested creep strain of the three-layered CGS (L1, L2 and L3) as well as its predicted ε-t curves with consideration of different interface-layer numbers: N = 1 (L1 + L1/L2 + L2 + L3 or L1 + L2 + L2/L3 + L3) and N = 2 (L1 + L1/L2 + L2 + L2/L3 + L3), where each interface-layer height is h = 0.3 H (H = 200 mm for the three-layered CGS, the total specimen height is 600 mm). Similarly, they all have the same three stages as those both the one-layered and two-layered CGS: instantaneous deformation, attenuation deformation, and stable creep deformation.

Comparing the Creep Curves of Two-Layered and Three-Layered CGS
The above triaxial creep test results of the one-layered, two-layered, and three-layered CGSs have proved the validity of increasing the degree of irregular pebble and height of interface-layer can decrease the creep deformation. For studying the effect of layer number on the creep curve of multilayered CGS, take two-layered CGS (L1 + L1/L2 + L2 or L2 + L2/L3 + L2) and three-layered (L1 + L1/L2 + L2 + L2/L1 + L1 or L2 + L2/L3 + L3 + L2/L3 + L2) CGS compacted by the same gradation of onelayered CGS (L1 and L2 or L2 and L3) under the same total height (600 mm) for comparison. Figure 14 shows the predicted creep curves of the two-layered CGS (each layer height H = 300 mm and interface-layer height h = 0.3 H = 90 mm) and three-layered CGS (each layer height H = 200 mm and interface-layer height h = 0.3 H = 60 mm). It is seen that the three-layered CGS has more interface-layer and smaller creep deformation than the two-layered CGS, since the existence of interface-layer decreases the creep deformation ( Figure 12). The more the interface-layer, the larger the total height of the interface-layer and thus the smaller the creep deformation. Obviously, the multi-layered CGS of more layer numbers can reduce settlement but also cause difficulty in field compaction and high cost. Both safety and economy must be considered in the design of a pavement or embankment.  The three-layered CGS has two interface-layers. The more the interface-layers are taken into account in calculation equation (Equation (19)), the predicted creep deformation becomes smaller, since it is known from Section 4.2 that the interface-layer with some height would decrease the deformation of multi-layered CGS, furthermore, three-layered CGS has two interface-layers with the height of 30% H, therefore, two interface-layers has the smallest deformation. The tested results agree well with the predicted results of the three-layered CGS with N = 1 or 2. This proves again the interface-layer really exists in the multi-layered CGS and the established interface-layer calculation model is reasonable and reliable. At least one interface-layer must be considered for predicting viscoelastic parameters of the multi-layered CGS.

Comparing the Creep Curves of Two-Layered and Three-Layered CGS
The above triaxial creep test results of the one-layered, two-layered, and three-layered CGSs have proved the validity of increasing the degree of irregular pebble and height of interface-layer can decrease the creep deformation. For studying the effect of layer number on the creep curve of multi-layered CGS, take two-layered CGS (L1 + L1/L2 + L2 or L2 + L2/L3 + L2) and three-layered (L1 + L1/L2 + L2 + L2/L1 + L1 or L2 + L2/L3 + L3 + L2/L3 + L2) CGS compacted by the same gradation of one-layered CGS (L1 and L2 or L2 and L3) under the same total height (600 mm) for comparison. Figure 14 shows the predicted creep curves of the two-layered CGS (each layer height H = 300 mm and interface-layer height h = 0.3 H = 90 mm) and three-layered CGS (each layer height H = 200 mm and interface-layer height h = 0.3 H = 60 mm). It is seen that the three-layered CGS has more interface-layer and smaller creep deformation than the two-layered CGS, since the existence of interface-layer decreases the creep deformation ( Figure 12). The more the interface-layer, the larger the total height of the interface-layer and thus the smaller the creep deformation. Obviously, the multi-layered CGS of more layer numbers can reduce settlement but also cause difficulty in field compaction and high cost. Both safety and economy must be considered in the design of a pavement or embankment.
interface-layer and smaller creep deformation than the two-layered CGS, since the existence of interface-layer decreases the creep deformation ( Figure 12). The more the interface-layer, the larger the total height of the interface-layer and thus the smaller the creep deformation. Obviously, the multi-layered CGS of more layer numbers can reduce settlement but also cause difficulty in field compaction and high cost. Both safety and economy must be considered in the design of a pavement or embankment.

Conclusions
1) The interface-layer of viscoelasticity and the actual shape of large-particle inclusion were firstly considered and a new interface-layer method was proposed to predict viscoelastic parameters of multi-layered CGS based on homogenization method and elastic-viscoelastic corresponding principle. The tested creep curves of the multi-layered CGS agreed well with the predicted ones and can prove the existence of the interface-layer and verify the validity of this new method. 2) For the one-layered CGS, the creep deformation decreased as the shape factor (ρ) of pebble

Conclusions
(1) The interface-layer of viscoelasticity and the actual shape of large-particle inclusion were firstly considered and a new interface-layer method was proposed to predict viscoelastic parameters of multi-layered CGS based on homogenization method and elastic-viscoelastic corresponding principle. The tested creep curves of the multi-layered CGS agreed well with the predicted ones and can prove the existence of the interface-layer and verify the validity of this new method. (2) For the one-layered CGS, the creep deformation decreased as the shape factor (ρ) of pebble inclusion increased and ρ was in the range of 1-1.8. (3) The viscoelastic parameters of the interface-layer were smaller than the average values of one-layered CGS which consisted of the interface-layer. (4) For the two-layered CGS, the creep deformation decreased with the increase of the interface-layer height (h) and the suitable interface-layer height was 20-30% as much as the height of one-layered CGS. (5) For the three-layered CGS, the creep deformation decreased with the increase of the interface-layer number (N) and at least one interface-layer must be taken into account. (6) For the safety point of view, it is better to use a high degree of irregular pebble and a large difference in gradation of CGS and more layer numbers to reduce the creep deformation of multi-layered coarse-grained soil. Construction costs should be taken into account in design of pavement or embankment.