Resilient Modulus Behavior and Prediction Models of Unbound Permeable Aggregate Base Materials Derived from Tunneling Rock Wastes

Tunneling rock wastes (TRWs), which are often open- or gap-graded in nature, have been increasingly recycled and reused for sustainable construction of unbound permeable aggregate base (UPAB) courses with high porosity and desired drainability. However, there is still a lack of sufficient understanding of long-term mechanical stability of such TRW materials subjected to repeated applications of moving wheel loads. This paper aimed to characterize and predict resilient modulus (Mr) behavior of the TRW materials used in unbound permeable aggregate base applications. To achieve this goal, five different UPAB gradations were designed based on the gravel-to-sand ratio (G/S) concept. In order to study their Mr behavior, the laboratory repeated load triaxial tests were conducted under different combinations of confining pressure and deviator stress as controlled by the levels of the shear stress ratio (SSR). The prediction accuracy of fourteen classical Mr prediction models was comparatively analyzed, from which the improved Mr prediction model incorporating gradation and stress variables was proposed for TRW-derived UPAB materials and further validated by external database accordingly. The results show that under the same G/S value and confining pressure level, the higher the SSR is, the greater the final Mr values are, and the more significant the effect of G/S on Mr is. Under the same SSR level, the increase of confining pressure alleviates the effect of G/S on Mr. There appears to exist an optimal G/S value of around 1.6–1.8 that yields the best Mr behavior of the TRW-derived UPAB materials studied. The improved Mr prediction model was verified extensively to be universally applicable. It can potentially contribute to balancing long-term mechanical stability and drainability of TRW-derived UPAB materials through gradation optimization. The findings could provide a theoretical basis and technical reference for cost-effective and sustainable applications of UPAB materials derived from TRWs.


Introduction
In recent years, the frequency of occurrence and severity of damage of catastrophic weather events, including heavy precipitation and urban flooding (known as the El Niño phenomenon), have increased significantly around the global [1][2][3]. In order to mitigate flooding risks and related damages, the initiative of sponge cities has been widely enforced, where the construction of permeable pavements has emerged as one of the effective alternatives to enhance drainage capacity of urban road networks [4][5][6]. As compared to traditional impervious or low-permeability pavement structures constructed with densegraded materials, permeable pavement layers constructed with open-graded materials

Materials Tested
The unbound aggregate materials used in the laboratory tests, as shown in Figure 1a, were sampled from a disposal facility of tunneling rock wastes in the suburb of Changsha, Hunan Province, and the lithology of these materials is limestone. A total of five different gradations were designed and controlled using the gravel-to-sand ratio (G/S) concept proposed previously [53] to represent UPAB materials with varying particle-packing structures and thus permeability levels. It is worth noting that the G/S concept was proposed based on the packing theory and has been verified extensively for its applicability to unbound granular materials, of which more details can be found elsewhere [53]. The G/S parameter calculated accordingly via Equations (1) and (2) was found to control key mechanical properties (e.g., shear strength and permanent deformation) of unbound aggregates reasonably well [53]. Note that the parameters D max and n in Equation (1) are fitted from the gradation data by the Talbot function in Equation (2), where all the sieve sizes and related percentages of materials passing those sieves ought to be used.
where P i is the percentage of materials passing the i-th sieve (%); d i is the opening aperture size of the i-th sieve (mm); D max is the maximum particle size (mm); and n is the shape parameter of the particular gradation curve, and it can be calculated from the given G/S value as per Equation (1), or be fitted from the gradation data by the Talbot function in Equation (2). It is worth noting that the n value was obtained by using the former method in this paper. In order to follow related Chinese design codes [54] and avoid size effects for triaxial tests [55,56], the maximum particle size was determined to be 26.5 mm. In this study, a suitable gradation band with a maximum particle size of 26.5 mm was chosen, which is also qualified for use as unbound permeable aggregate gradations according to the Chinese technical specifications for highway pavement base construction [54,57]. Within this chosen gradation band, the corresponding G/S values range from 1.57-2.8. As plotted in Figure 1b, a total of five different gradation curves were selected from this gradation band to prepare TRW-derived UPAB specimens, and their G/S values were calculated to be 1.0, 1.6, 1.8, 2.0, and 2.5, respectively. The gradation curves with the G/S values of 1.6, 1.8, 2.0, and 2.5 represented low-, medium-, and high-level permeability, whereas the one with the G/S value of 1.0 served as the baseline against which the other four were compared due to its extremely low permeability. The G/S range of 1.0-2.5 studied is believed to well cover typical gradations of unbound aggregates used in pavement base construction practices. The related parameters of the gradation design scheme are summarized in Table 1. Note: C c is the coefficient of curvature; C u is the coefficient of uniformity.

Laboratory Compaction Tests
The TRW-derived unbound aggregate specimens with the G/S values of 1.0, 1.6, 1.8, 2.0, and 2.5 were prepared at five different initial moisture content levels of 3%, 4%, 4.5%, 5%, and 6% to carry out laboratory Proctor compaction tests, respectively. The inner diameter and height of the compaction mold were 150 mm and 120 mm, respectively. According to the Chinese geotechnical standards [57], the specimens were compacted in 3 sub-layers, with each sub-layer subjected to 98 blows. The applied compaction energy per unit volume was 2677 kJ/m 3 , which is similar to the 2693.3 kJ/m 3 of the modified Proctor compaction according to AASHTO T180 [58]. The compaction curves of gravitational moisture content versus achieved dry density for specimens with different G/S values were obtained and are plotted in Figure 2, from which their optimal moisture content and maximum dry density values were determined accordingly.

Laboratory Compaction Tests
The TRW-derived unbound aggregate specimens with the G/S values of 1.0, 1.6, 1.8, 2.0, and 2.5 were prepared at five different initial moisture content levels of 3%, 4%, 4.5%, 5%, and 6% to carry out laboratory Proctor compaction tests, respectively. The inner diameter and height of the compaction mold were 150 mm and 120 mm, respectively. According to the Chinese geotechnical standards [57], the specimens were compacted in 3 sub-layers, with each sub-layer subjected to 98 blows. The applied compaction energy per unit volume was 2677 kJ/m 3 , which is similar to the 2693.3 kJ/m 3 of the modified Proctor compaction according to AASHTO T180 [58]. The compaction curves of gravitational moisture content versus achieved dry density for specimens with different G/S values were obtained and are plotted in Figure 2, from which their optimal moisture content and maximum dry density values were determined accordingly.
It can be seen from Figure 2 that the compaction curves of the unbound aggregate specimens with different G/S values all exhibit similar trends, i.e., the achieved dry density first increases and then decreases with increasing moisture content, with clear peaks observed. The optimum moisture content value of the specimens with different G/S values is around 4.5%, according to Figure 2. However, the maximum dry density values of specimens with different G/S values are also different. As the G/S value gradually increases from 1.0 to 2.5, the maximum dry density of the specimens first increases and then decreases. There seems to exist a peak value of the maximum dry density of the specimens, of which the corresponding G/S value is between 1.6 and 1.8. This is consistent with the result reported by Zhou et al. [59]. This is expected, as fine fractions dominate for low G/S values and coarse fractions dominate for high G/S values. Therefore, specimens with either low or high G/S values are relatively more difficult to compact, thus justifying the existence of the peak value of maximum dry density. It can be seen from Figure 2 that the compaction curves of the unbound aggregate specimens with different G/S values all exhibit similar trends, i.e., the achieved dry density first increases and then decreases with increasing moisture content, with clear peaks observed. The optimum moisture content value of the specimens with different G/S values is around 4.5%, according to Figure 2. However, the maximum dry density values of specimens with different G/S values are also different. As the G/S value gradually increases from 1.0 to 2.5, the maximum dry density of the specimens first increases and then decreases. There seems to exist a peak value of the maximum dry density of the specimens, of which the corresponding G/S value is between 1.6 and 1.8. This is consistent with the result reported by Zhou et al. [59]. This is expected, as fine fractions dominate for low G/S values and coarse fractions dominate for high G/S values. Therefore, specimens with either low or high G/S values are relatively more difficult to compact, thus justifying the existence of the peak value of maximum dry density.

Monotonic Triaxial Compression Tests
The triaxial test equipment is of the model No. STD and manufactured by the XIAN LETRY Company, Xi'an, China. The dimensions of the laboratory triaxial specimens were 100 mm in diameter and 200 mm in height with the diameter-to-height ratio of 1:2. The unbound aggregate specimens with different G/S values were fabricated at their respective optimal conditions (defined by the maximum dry density and related optimal moisture content) and then compacted in 5 sub-layers to achieve the degree of compaction of at least 98%. The degree of compaction is defined as the ratio of achieved dry density to the related maximum dry density at the same or similar compaction energy level. The

Monotonic Triaxial Compression Tests
The triaxial test equipment is of the model No. STD and manufactured by the XIAN LETRY Company, Xi'an, China. The dimensions of the laboratory triaxial specimens were 100 mm in diameter and 200 mm in height with the diameter-to-height ratio of 1:2. The unbound aggregate specimens with different G/S values were fabricated at their respective optimal conditions (defined by the maximum dry density and related optimal moisture content) and then compacted in 5 sub-layers to achieve the degree of compaction of at least 98%. The degree of compaction is defined as the ratio of achieved dry density to the related maximum dry density at the same or similar compaction energy level. The compacted specimens were restrained with 0.5 mm thick latex membranes and then transferred to the triaxial chamber for subsequent laboratory monotonic and repeated load triaxial compression tests.
The triaxial testing apparatus is shown in Figure 3. It is capable of applying both monotonic and repeated axial loads of different amplitude and/or frequency levels. The vertical deformation is recorded externally with the linear variable differential transducers (LVDTs). Each aggregate specimen was loaded under three different levels of confining pressure (σ 3 ) of 50 kPa, 100 kPa, and 150 kPa, respectively. The half-sine waveform with a loading frequency of 5 Hz was adopted for deviator stress (σ d ), as shown in Figure 4. The shear strength parameters (i.e., the apparent cohesion c and the angle of internal friction φ) of the specimens were determined by fitting the results of monotonic load triaxial compression tests against the classical Mohr-Coulomb failure criterion.

Repeated Load Triaxial Tests
The repeated load triaxial (RLT) tests were further conducted to study the resilient and plastic strain behavior of TRW-derived UPAB materials, i.e., resilient modulus and permanent deformation characteristics. The shear stress ratio (SSR), of which the definition is illustrated in Figure 5, was used to control the amplitudes of cyclic deviator stress (σd). The related calculation formulas are shown in Equations (3)-(6) sequentially [23]. To study the effect of stress level on resilient modulus of UPAB materials, four different SSR levels were designed for repeated load triaxial tests, i.e., 0.3, 0.5, 0.7, and 0.9, respectively. Therefore, each of the five different gradations corresponds to a total of 12 different combinations of stress states (i.e., 3 different σ3 levels times 4 different SSR levels), as listed in Table 2. The permanent deformation of the UPAB materials tended to be stable after 2000 load applications [53,55]. The dynamic properties of the UPAB materials were basically stable when the number of load applications exceeded 5000 [56]. By considering the stress levels applied to the specimens and the capacity of the triaxial testing apparatus, it was determined in this study that the maximum number of repeated load applications was set as 50,000. The specimens of five different gradations subjected to 50,000 load applications

Repeated Load Triaxial Tests
The repeated load triaxial (RLT) tests were further conducted to study the resilient and plastic strain behavior of TRW-derived UPAB materials, i.e., resilient modulus and permanent deformation characteristics. The shear stress ratio (SSR), of which the definition is illustrated in Figure 5, was used to control the amplitudes of cyclic deviator stress (σd). The related calculation formulas are shown in Equations (3)-(6) sequentially [23]. To study the effect of stress level on resilient modulus of UPAB materials, four different SSR levels were designed for repeated load triaxial tests, i.e., 0.3, 0.5, 0.7, and 0.9, respectively. Therefore, each of the five different gradations corresponds to a total of 12 different combinations of stress states (i.e., 3 different σ3 levels times 4 different SSR levels), as listed in Table 2. The permanent deformation of the UPAB materials tended to be stable after 2000 load applications [53,55]. The dynamic properties of the UPAB materials were basically stable when the number of load applications exceeded 5000 [56]. By considering the stress levels applied to the specimens and the capacity of the triaxial testing apparatus, it was determined in this study that the maximum number of repeated load applications was set as 50,000. The specimens of five different gradations subjected to 50,000 load applications

Repeated Load Triaxial Tests
The repeated load triaxial (RLT) tests were further conducted to study the resilient and plastic strain behavior of TRW-derived UPAB materials, i.e., resilient modulus and permanent deformation characteristics. The shear stress ratio (SSR), of which the definition is illustrated in Figure 5, was used to control the amplitudes of cyclic deviator stress (σ d ). The related calculation formulas are shown in Equations (3)-(6) sequentially [23]. To study the effect of stress level on resilient modulus of UPAB materials, four different SSR levels were designed for repeated load triaxial tests, i.e., 0.3, 0.5, 0.7, and 0.9, respectively. Therefore, each of the five different gradations corresponds to a total of 12 different combinations of stress states (i.e., 3 different σ 3 levels times 4 different SSR levels), as listed in Table 2. The permanent deformation of the UPAB materials tended to be stable after 2000 load applications [53,55]. The dynamic properties of the UPAB materials were basically stable when the number of load applications exceeded 5000 [56]. By considering the stress levels applied to the specimens and the capacity of the triaxial testing apparatus, it was determined in this study that the maximum number of repeated load applications was set as 50,000. The specimens of five different gradations subjected to 50,000 load applications under 50 kPa confining pressure and 93. -kPa deviator stress were photographed and shown as examples in Figure 6.
where M r is resilient modulus, ε r is resilient (or recoverable) axial strain, ε p is irreversible (or plastic) axial strain, λ is damping ratio, S cycle is the enclosed area of the axial stress-axial strain hysteresis loop, and S OAA is the enclosed area of the triangle OAA , as shown in Figure 7.
where r M is resilient modulus, r ε is resilient (or recoverable) axial strain, p ε is irreversible (or plastic) axial strain, λ is damping ratio, is the enclosed area of the axial stress-axial strain hysteresis loop, and OAA S ′ is the enclosed area of the triangle OAA′ , as shown in Figure 7.       values by using the obtained shear strength parameters (c and ϕ).    Figure 8 shows the variation of accumulated axial plastic strain with the number of load applications (N) for TRW-derived UPAB specimens with different G/S levels under 50 kPa confining pressure. It is worth noting that the amplitudes of deviator stress (σ d ) applied on the specimens with different G/S values are different even under the same combination of confining pressure and SSR due to their differences in shear strength. The accumulated axial plastic strain of the specimens increases rapidly with increasing number of load applications (N), while the rate of increase decreases during the initial loading stage. The accumulated axial plastic strain tends to be stable approximately at lower SSR levels (≤0.5), while it accumulates rapidly with increasing number of load applications (N) at higher SSR levels (≥0.7). The accumulated axial plastic strain increases significantly with increasing SSR under the same confining pressure and gradation. Figure 8 shows the variation of accumulated axial plastic strain with the number of load applications (N) for TRW-derived UPAB specimens with different G/S levels under 50 kPa confining pressure. It is worth noting that the amplitudes of deviator stress (σd) applied on the specimens with different G/S values are different even under the same combination of confining pressure and SSR due to their differences in shear strength. The accumulated axial plastic strain of the specimens increases rapidly with increasing number of load applications (N), while the rate of increase decreases during the initial loading stage. The accumulated axial plastic strain tends to be stable approximately at lower SSR levels (≤0.5), while it accumulates rapidly with increasing number of load applications (N) at higher SSR levels (≥0.7). The accumulated axial plastic strain increases significantly with increasing SSR under the same confining pressure and gradation. In order to further study the influence of the gradation parameter G/S on the accumulated axial plastic strain (ε1p acc ) of TRW-derived UPAB materials under different combinations of σ3 and SSR, the ε1p acc values calculated upon the completion of laboratory repeated load triaxial (RLT) tests (i.e., N = 50,000) are plotted for comparison in Figure 9. It can be seen from Figure 9 that under the same combination of σ3 and SSR, the ultimate accumulated axial plastic strain first increases and then decreases with increasing G/S level. This implies that the ultimate accumulated axial plastic strain reaches its minimum and maximum values at G/S = 1.6 and G/S = 2.5, respectively. Therefore, the trend of the In order to further study the influence of the gradation parameter G/S on the accumulated axial plastic strain (ε 1p acc ) of TRW-derived UPAB materials under different combinations of σ 3 and SSR, the ε 1p acc values calculated upon the completion of laboratory repeated load triaxial (RLT) tests (i.e., N = 50,000) are plotted for comparison in Figure 9. It can be seen from Figure 9 that under the same combination of σ 3 and SSR, the ultimate accumulated axial plastic strain first increases and then decreases with increasing G/S level. This implies that the ultimate accumulated axial plastic strain reaches its minimum and maximum values at G/S = 1.6 and G/S = 2.5, respectively. Therefore, the trend of the ultimate accumulated axial plastic strain varying with G/S is consistent with previously presented trends of maximum dry density and shear strength varying with G/S, respectively. Figure 10 shows the variations of M r with the number of load applications (N) for TRW-derived UPAB specimens with different G/S levels under 50 kPa confining pressure. It can be seen from Figure 10 that M r increases rapidly with increasing number of load applications (N) at the initial loading stage until the gradually decreasing rate of increase becomes stable approximately at N = 1500. In addition, the M r values of UPAB specimens all increase significantly with increasing SSR level, but still vary among different G/S levels even under the same stress state. The UPAB specimen with G/S = 2.5 exhibits the lowest M r values, which may be attributable to its high fraction of coarse particles and large pores; in contrast, the UPAB specimen with G/S = 1.6 exhibits the greatest M r values due to its relatively densest packing structure. Figure 10 also shows that the higher the SSR level is, the more significant the influence of gradation on M r of UPAB specimens is, and the more pronounced the increases in M r become.

Resilient Modulus Characteristics
Materials 2022, 15, x FOR PEER REVIEW 11 ultimate accumulated axial plastic strain varying with G/S is consistent with previo presented trends of maximum dry density and shear strength varying with G/S, re tively.  Figure 10 shows the variations of Mr with the number of load applications (N TRW-derived UPAB specimens with different G/S levels under 50 kPa confining pres It can be seen from Figure 10 that Mr increases rapidly with increasing number of applications (N) at the initial loading stage until the gradually decreasing rate of inc becomes stable approximately at N = 1500. In addition, the Mr values of UPAB specim all increase significantly with increasing SSR level, but still vary among different G/S els even under the same stress state. The UPAB specimen with G/S = 2.5 exhibits the lo Mr values, which may be attributable to its high fraction of coarse particles and pores; in contrast, the UPAB specimen with G/S = 1.6 exhibits the greatest Mr values to its relatively densest packing structure. Figure 10 also shows that the higher the level is, the more significant the influence of gradation on Mr of UPAB specimens is the more pronounced the increases in Mr become.   Figure 10 shows the variations of Mr with the number of load applications (N TRW-derived UPAB specimens with different G/S levels under 50 kPa confining press It can be seen from Figure 10 that Mr increases rapidly with increasing number of applications (N) at the initial loading stage until the gradually decreasing rate of incr becomes stable approximately at N = 1500. In addition, the Mr values of UPAB specim all increase significantly with increasing SSR level, but still vary among different G/S els even under the same stress state. The UPAB specimen with G/S = 2.5 exhibits the lo Mr values, which may be attributable to its high fraction of coarse particles and l pores; in contrast, the UPAB specimen with G/S = 1.6 exhibits the greatest Mr values to its relatively densest packing structure. Figure 10 also shows that the higher the level is, the more significant the influence of gradation on Mr of UPAB specimens is, the more pronounced the increases in Mr become.   Figure 11 that under the same combination of σ 3 and SSR, M r values first increase and then decrease with increasing G/S level, thus clearly indicating that the highest M r is reached at G/S = 1.6 and the lowest M r at G/S = 2.5. Therefore, it appears that the optimal G/S value of 1.6-1.8 exists where the maximum resilient modulus (or the best resilient behavior) can be expected for UPAB materials. Under the same confining pressure level of 150 kPa and the varying SSR levels of 0.3, 0.5, 0.7, and 0.9, M r values of the specimen with G/S = 1.6 increase by 4.53%, 4.99%, 5.95%, and 7.18% as compared to their counterparts of the specimens with G/S = 2.5, respectively. The influence of gradation on M r characteristics of UPAB materials appears more pronounced for higher SSR levels. In addition, under the same SSR level of 0.5, as confining pressure increases, M r values of the specimens with G/S = 1.6 increase by 7.64%, 5.24%, and 4.99% as compared to their counterparts of the specimens with G/S = 2.5, respectively. This implies that under the same SSR level, increasing confining pressure would alleviate the effect of gradation on resilient modulus of unbound permeable aggregate base materials, which is expected by the common engineering sense.   Figure 11 that under the same combination of σ3 and SSR, Mr values first increase and then decrease with increasing G/S level, thus clearly indicating that the highest Mr is reached at G/S = 1.6 and the lowest Mr at G/S = 2.5. Therefore, it appears that the optimal G/S value of 1.6-1.8 exists where the maximum resilient modulus (or the best resilient behavior) can be expected for UPAB materials. Under the same confining pressure level of 150 kPa and the varying SSR levels of 0.3, 0.5, 0.7, and 0.9, Mr values of the specimen with G/S = 1.6 increase by 4.53%, 4.99%, 5.95%, and 7.18% as compared to their counterparts of the specimens with G/S = 2.5, respectively. The influence of gradation on Mr characteristics of UPAB materials appears more pronounced for higher SSR levels. In addition, under the same SSR level of 0.5, as confining pressure increases, Mr values of the specimens with G/S = 1.6 increase by 7.64%, 5.24%, and 4.99% as compared to their counterparts of the specimens with G/S = 2.5, respectively. This implies that under the same SSR level, increasing confining pressure would alleviate the effect of gradation on resilient modulus of unbound permeable aggregate base materials, which is expected by the common engineering sense.

Comparative Assessment of Different Calculation Methods of Representative Mr Values
In the past few decades, resilient modulus prediction models for subgrade soils and unbound granular materials have evolved dramatically from primitive to sophisticated ones with key physical properties and the shearing and confinement effects considered sequentially. Among the existing Mr models in the literature, fourteen classical ones were shortlisted for comparative assessment in this study, as listed in Table 3. Since the matrix suction was not scoped into or measured during the laboratory testing program, those Mr prediction models incorporating matrix-suction-related terms were excluded from the scope of this study.

Comparative Assessment of Different Calculation Methods of Representative M r Values
In the past few decades, resilient modulus prediction models for subgrade soils and unbound granular materials have evolved dramatically from primitive to sophisticated ones with key physical properties and the shearing and confinement effects considered sequentially. Among the existing M r models in the literature, fourteen classical ones were shortlisted for comparative assessment in this study, as listed in Table 3. Since the matrix suction was not scoped into or measured during the laboratory testing program, those M r prediction models incorporating matrix-suction-related terms were excluded from the scope of this study.
In order to study the influence of different calculation methods of the representative M r values on the prediction accuracy of the M r models, three different methods for calculating the representative M r values were comparatively assessed. The first method averaged M r values from the last 5 of the first 1500 load applications. The second method averaged M r values from the last 10 of the first 10,000 load applications, whereas the third one averaged M r values from the last 100 of the first 50,000 load applications. The entire sets of such representative M r data calculated by using three different methods for UPAB specimens tested under different combinations of G/S levels and stress states were fitted against each of the M r prediction models listed in Table 3. The resulting goodness-of-fit results (e.g., R 2 values) obtained for the fourteen different M r prediction models according to the three different calculation methods of the representative Mr values are summarized in Table 3 as well. Notes: M r denotes resilient modulus (MPa); k 1 , k 2 , k 3 , k 4 , and k 5 denote model coefficients to be fitted; θ denotes the first principal stress invariant (kPa), i.e., θ = σ 1 + σ 2 + σ 3 ; τ oct denotes the octahedral shear stress, i.e., 3 σ d in triaxial stress space where σ 2 = σ 3 ); and P a denotes the reference stress that is generally taken as the standard atmospheric pressure of 100 kPa. Table 3 shows that the R 2 values of the fourteen different M r models for the calculation method I are slightly greater than those for the calculation methods II and III. Overall speaking, little difference is observed among such three different methods in terms of the prediction accuracy, thus indicating that any one of them is suitable for calculating representative M r values in this study. By comparing the R 2 values of the fourteen different M r models, it further proves that the addition of both standard atmospheric pressure (Pa) and the constant of 1 has little effect on the improvement of prediction accuracy, and that the inclusion of additional model coefficients k 4 and k 5 is not effective for improving the prediction accuracy of the fourteen different M r prediction models.

Comparative Assessment of Existing Classical M r Prediction Models
The M r prediction models listed in Table 3 were employed to fit the laboratorymeasured representative M r results under different combinations of G/S levels and stress states. The goodness-of-fit indicators including the coefficient of determination (R 2 ), ad-justed R 2 (adj. R 2 ), and root mean square error (RMSE) were calculated for comparative assessment of such different M r models. Since the axial strain (or deformation) of the specimens recorded during each load application of laboratory RLT tests consists of both resilient and plastic parts, the resilient strain induced by each load application can be used to calculate M r and evaluate its variation with the number of load applications under different combinations of G/S and stress levels, whereas the plastic strain can be used to evaluate permanent deformation accumulation behavior. To facilitate such comparisons, the above-mentioned method I was adopted herein to calculate the representative M r values of each specimen. The results of such representative M r values were then fitted by each of the M r prediction models listed in Table 3 to assess their applicability and accuracy of prediction. The specimens with G/S = 1.0 were taken as examples to illustrate this process and the related results in the following subsections. Although not presented herein for brevity, the comparison results of laboratory-measured versus model-predicted representative M r values for UPAB specimens with other G/S values than 1.0 (i.e., G/S = 1.6, 1.8, 2.0, or 2.5) also exhibit similar trends and observations to be detailed subsequently. Figure 12 plots the 12 sets of laboratory-measured representative M r data for the specimens with G/S = 1.0, along with the results predicted by classical univariate M r models. It can be seen from Figure 12 that the R 2 value is only 0.464 for both Model No. 1 and Model No. 2 incorporating the first principal stress invariant (θ) only. This indicates that the prediction accuracy of both Models No. 1 and No. 2 is poor, and that normalizing the stress variable by the standard atmospheric pressure barely improves model accuracy. In contrast, the R 2 value of Model No. 3 incorporating deviatoric stress (σ d ) only reaches 0.962, thus indicating that the prediction accuracy is relatively satisfactory. It further proves that the correlation between M R and σ d is more significant than θ for UPAB materials studied. inclusion of additional model coefficients k4 and k5 is not effective for improving the prediction accuracy of the fourteen different Mr prediction models.

Comparative Assessment of Existing Classical Mr Prediction Models
The Mr prediction models listed in Table 3 were employed to fit the laboratory-measured representative Mr results under different combinations of G/S levels and stress states. The goodness-of-fit indicators including the coefficient of determination (R 2 ), adjusted R 2 (adj. R 2 ), and root mean square error (RMSE) were calculated for comparative assessment of such different Mr models. Since the axial strain (or deformation) of the specimens recorded during each load application of laboratory RLT tests consists of both resilient and plastic parts, the resilient strain induced by each load application can be used to calculate Mr and evaluate its variation with the number of load applications under different combinations of G/S and stress levels, whereas the plastic strain can be used to evaluate permanent deformation accumulation behavior. To facilitate such comparisons, the abovementioned method I was adopted herein to calculate the representative Mr values of each specimen. The results of such representative Mr values were then fitted by each of the Mr prediction models listed in Table 3 to assess their applicability and accuracy of prediction. The specimens with G/S = 1.0 were taken as examples to illustrate this process and the related results in the following subsections. Although not presented herein for brevity, the comparison results of laboratory-measured versus model-predicted representative Mr values for UPAB specimens with other G/S values than 1.0 (i.e., G/S = 1.6, 1.8, 2.0, or 2.5) also exhibit similar trends and observations to be detailed subsequently. Figure 12 plots the 12 sets of laboratory-measured representative Mr data for the specimens with G/S = 1.0, along with the results predicted by classical univariate Mr models. It can be seen from Figure 12 that the R 2 value is only 0.464 for both Model No. 1 and Model No. 2 incorporating the first principal stress invariant (θ) only. This indicates that the prediction accuracy of both Models No. 1 and No. 2 is poor, and that normalizing the stress variable by the standard atmospheric pressure barely improves model accuracy. In contrast, the R 2 value of Model No. 3 incorporating deviatoric stress ( ) only reaches 0.962, thus indicating that the prediction accuracy is relatively satisfactory. It further proves that the correlation between MR and is more significant than for UPAB materials studied. The comparisons of laboratory-measured representative M r values against predicted values by different bivariate (i.e., σ 3 and σ d ) M r prediction models are shown in Figure 13 for UPAB specimens with G/S = 1.0. It can be observed from Figure 13  13 for UPAB specimens with G/S = 1.0. It can be observed from Figure 13 that the R 2 values of those bivariate Mr prediction models (i.e., Model No. 4,Model No. 5,and Model No. 6) incorporating both and are all above 0.95 and insignificantly different from each other. Due to the dimension inconsistency of Model No. 4, it was abandoned and excluded from further analysis. By comparing the R 2 values between Mr prediction models No. 5 and No. 6, it can be found that the addition of the constant of 1 to stress variables in Model No. 6 results in almost negligible difference in prediction accuracy.  Similarly, the laboratory-measured representative Mr values and those predicted by other prediction models incorporating the first principal stress invariant (θ) and the octahedral shear stress ( ) (i.e., Models No. 10 to No. 14) are compared and plotted in Figure  15 for the UPAB specimens with G/S = 1.0. It can be seen from Figure 15 that the R 2 values of Models No. 10,No. 11,and No. 12 are consistent, thus indicating that the addition of standard atmospheric pressure (Pa) and the constant of 1 has little effect on improving the prediction accuracy. However, the R 2 values of Models No. 14 and No. 15 slightly decrease as compared to those of Models No. 10 to No. 12, indicating that the addition of model coefficients k4 and k5 is not effective for accuracy improvement.   Similarly, the laboratory-measured representative M r values and those predicted by other prediction models incorporating the first principal stress invariant (θ) and the octahedral shear stress (τ oct ) (i.e., Models No. 10 to No. 14) are compared and plotted in Figure 15 for the UPAB specimens with G/S = 1.0. It can be seen from Figure 15 that the R 2 values of Models No. 10,No. 11,and No. 12 are consistent, thus indicating that the addition of standard atmospheric pressure (Pa) and the constant of 1 has little effect on improving the prediction accuracy. However, the R 2 values of Models No. 14 and No. 15 slightly decrease as compared to those of Models No. 10 to No. 12, indicating that the addition of model coefficients k 4 and k 5 is not effective for accuracy improvement.

Selection of the Optimal Classical M r Prediction Model
The aforementioned M r prediction models with R 2 values greater than 0.95 were further shortlisted to perform additional assessment. Due to the dimension inconsistency existing in the formulas of Models No. 4 and No. 7, they were excluded from further analysis. As presented previously, the addition of standard atmospheric pressure (Pa) and the constant of 1 has little effect on the improvement of prediction accuracy, and the inclusion of additional model coefficients k 4 and k 5 reduces the prediction accuracy. Therefore, the focus was to further compare Models No. 3, No. 5, No. 8, and No. 10 through comprehensive statistical analysis. Table 4 summarizes the statistics of the residual eigenvalues between the laboratory-measured representative M r values and their counterparts predicted by Models No. 3, No. 5, No. 8, and No. 10 for the entire UPAB specimens, respectively. Model No. 8 is observed to have the lowest residual sum of squares, variance, root mean square error (RMSE), and Akaike information criterion (AIC) index, as compared to the other three models. In addition, the skewness of Model No. 8 is the lowest and its kurtosis is relatively large, thus indicating that its residual distribution is the closest to the normal distribution and its goodness-of-fit is the best. Therefore, in terms of the overall prediction accuracy, Model No. 8 is the best among the four M r models evaluated for the UPAB materials studied. Figure 16 plots the P-P diagrams of the normalized residuals between the measured and predicted representative M r values. Note that the P-P diagrams were drawn according to the relations between the cumulative proportions of predicted M r values and those of the normal distributions. That is, when the dataset conforms to the normal distribution, each data point approximately falls on the straight line of normal distribution in the P-P graph. It can be found from Figure 16 that the normal P-P diagrams of the standardized residuals generated by the four M r prediction models all coincide approximately with the straight line. The P-P diagram of Model No. 8 is much closer to the straight line than those of other models. This further confirms that the goodness-of-fit of Model No. 8 is the best among those four models evaluated.  Hence, the optimal prediction model was determined to be Model No. 8. It was then used to fit the measured representative Mr data of the UPAB specimens tested under different combinations of G/S levels and stress states in the laboratory repeated load triaxial tests. The model coefficients k1, k2, and k3 and related adjusted coefficient of determination (Adj. R 2 ) were determined for each G/S level accordingly and are listed in Table 5. It can be seen from Table 5 that as the G/S value increases, the Adj. Hence, the optimal prediction model was determined to be Model No. 8. It was then used to fit the measured representative M r data of the UPAB specimens tested under different combinations of G/S levels and stress states in the laboratory repeated load triaxial tests. The model coefficients k 1 , k 2, and k 3 and related adjusted coefficient of determination (Adj. R 2 ) were determined for each G/S level accordingly and are listed in Table 5. It can be seen from Table 5 that as the G/S value increases, the Adj.

Modification of the Optimal M r Prediction Model
As shown previously, the values of the three model coefficients k 1 , k 2, and k 3 of the optimal M r prediction Model No. 8 are considerably different from each other. In order to further address the influences of gradation variation (quantified by varying G/S levels), statistical correlation analysis was performed between the gradation control parameter G/S and the fitted model coefficients k 1 , k 2 , and k 3 . Different types of fitting functions (e.g., linear, quadratic, power, and logarithmic) were attempted, respectively. It was found by trial and error that there exists a statistically significant quadratic function relating model coefficient k 1 and the gradation parameter G/S, and a power function relating k 2 and G/S. However, k 3 and G/S are not statistically correlated. Therefore, Equations (8) and (9) were obtained accordingly, thus leading to the modified version of the M r prediction Model No. 8, as formulated in Equation (10). Note that the model coefficient k 3 in the M r prediction Model No. 8 varied insignificantly among different G/S levels and thus was eventually set as the constant of 0.405, the average k 3 value from different G/S levels.
0.405 (10) where M r is resilient modulus, k 1 and k 2 are model coefficients to be regressed, G/S is the gradation control parameter, σ d is deviatoric stress (kPa), θ is the first principal stress invariant (kPa), and Pa is the standard atmospheric pressure. In order to substantiate the validity and accuracy of the above-modified optimal M r prediction model (see Equation (10)), its predicted representative M r values were compared against the laboratory-measured counterparts of a total of 60 UPAB specimens, as shown in Figure 17b. It can be clearly seen from Figure 17b that the prediction errors increase with increasing M r values, and that the prediction errors for greater M r values are mainly attributable to the specimens with G/S values of 2 and 2.5. Therefore, it becomes imperative to further improve Equation (10) by dividing it by the arithmetric square root of the gradation parameter G/S, as shown in Equation (12).
The representative M r values of all the UPAB specimens predicted by Equation (11) were compared against their laboratory-measured counterparts accordingly. Figure 17c shows that the resulting R 2 value of Equation (11) increases from 0.73 to 0.83 as compared to Equation (10), indicating that the modification made to Equation (10) (i.e., Equation (11)) is effective in terms of enhancing the M r prediction accuracy for UPAB materials. The prediction accuracy of the original M r Model No. 8 is shown in Figure 17a as well, where the R 2 value only reaches 0.53 and is much lower than that yielded by the modified model (i.e., Equation (11)). Therefore, the modified model formulated in Equation (11) was finalized as the optimal M r prediction model.
The representative Mr values of all the UPAB specimens predicted by Equation (11) were compared against their laboratory-measured counterparts accordingly. Figure 17c shows that the resulting R 2 value of Equation (11) increases from 0.73 to 0.83 as compared to Equation (10), indicating that the modification made to Equation (10) (i.e., Equation (11)) is effective in terms of enhancing the Mr prediction accuracy for UPAB materials. The prediction accuracy of the original Mr Model No. 8 is shown in Figure 17a as well, where the R 2 value only reaches 0.53 and is much lower than that yielded by the modified model (i.e., Equation (11)). Therefore, the modified model formulated in Equation (11) was finalized as the optimal Mr prediction model.

Parametric Sensitivity Analysis of the Final Optimal Mr Predictive Model
The sensitivity of Mr to gradation variations under different stress states was studied by using the final optimal Mr prediction model (i.e., Equation (11)). The results of such

Parametric Sensitivity Analysis of the Final Optimal M r Predictive Model
The sensitivity of M r to gradation variations under different stress states was studied by using the final optimal Mr prediction model (i.e., Equation (11)). The results of such parametric sensitivity analysis are plotted in Figure 18. It can be seen from Figure 18 that for G/S values within the range of 1.0-2.5, the representative M r value predicted by Equation (11) decreases with increasing G/S value under the same combination of confining pressure and deviatoric stress, that is, as the content of fine particles decreases, the representative M r values decrease. In addition, the gradation variation has a much greater influence on the representative M r for lower G/S values. This implies that the representative M r is more sensitive to the gradation variation for UPAB materials with higher fraction of fine particles. In contrast, the gradation variation has little effect on the representative M r for relatively large G/S value, i.e., the representative M r is less sensitive to the gradation variation for the UPAB materials with higher fraction of coarse particles. parametric sensitivity analysis are plotted in Figure 18. It can be seen from Figure 18 that for G/S values within the range of 1.0-2.5, the representative Mr value predicted by Equation (11) decreases with increasing G/S value under the same combination of confining pressure and deviatoric stress, that is, as the content of fine particles decreases, the representative Mr values decrease. In addition, the gradation variation has a much greater influence on the representative Mr for lower G/S values. This implies that the representative Mr is more sensitive to the gradation variation for UPAB materials with higher fraction of fine particles. In contrast, the gradation variation has little effect on the representative Mr for relatively large G/S value, i.e., the representative Mr is less sensitive to the gradation variation for the UPAB materials with higher fraction of coarse particles.

Verification of the Developed Resilient Modulus Prediction Model
The additional external resilient modulus database was collected from the literature [67,68] to further verify the previously developed and validated Mr prediction model (i.e., Equation (11)). The three different types of commonly-used unbound aggregate materials were tested, i.e., uncrushed gravel, crushed limestone, and crushed dolomite. Four different levels of fines content (i.e., materials passing No. 200 or 0.075 mm sieve) were designed for each type of unbound aggregate materials, i.e., 4%, 8%, 12%, and 18%. They correspond to the G/S values of 1.39, 1.43, 1.49, and 1.55, respectively. In addition, two different types of fines, i.e., plastic (clay with plasticity index of 10-20%) and nonplastic (mineral

Verification of the Developed Resilient Modulus Prediction Model
The additional external resilient modulus database was collected from the literature [67,68] to further verify the previously developed and validated M r prediction model (i.e., Equation (11)). The three different types of commonly-used unbound aggregate materials were tested, i.e., uncrushed gravel, crushed limestone, and crushed dolomite. Four different levels of fines content (i.e., materials passing No. 200 or 0.075 mm sieve) were designed for each type of unbound aggregate materials, i.e., 4%, 8%, 12%, and 18%. They correspond to the G/S values of 1.39, 1.43, 1.49, and 1.55, respectively. In addition, two different types of fines, i.e., plastic (clay with plasticity index of 10-20%) and nonplastic (mineral powder with plasticity index of near zero), were used for each type of unbound aggregate materials as well. Although each type of unbound aggregate materials was tested under three different levels of moisture content (i.e., 90% × w opt , w opt , and 110% × w opt ) to study the influence of moisture content, only the M r data measured at the optimal conditions were selected for use in the model verification, as the final optimal M r prediction model was previously developed using the M r data measured at the optimal conditions only in this study. Table 6 summarizes the basic physical properties of the three different types of unbound aggregate materials tested in the literature [67,68]. Notes: L denotes limestone; D denotes dolomite; G denotes gravel; P denotes plastic (clay); NP denotes nonplastic (mineral power); and the number in the symbol of each material denotes the fines content (i.e., percentage of materials passing No. 200 or 0.075 mm sieve), for example, "L4NP" denotes unbound crushed limestone aggregate materials with the nonplastic fines content of 4%.
By using the above-mentioned external M r datasets, Table 7 compares the laboratorymeasured M r values against those predicted by the final optimal model established in Equation (11) for each of the three different types of unbound aggregate materials. It can be seen from Table 7 that the M r prediction model established in this study for UPAB materials is also reasonably accurate for the external unbound aggregates with different combinations of mineralogy, fines content, and plasticity index of fines. Overall, it works better for either limestone aggregate materials or aggregate materials with plastic fines as compared to other materials. Note that the model coefficients in Equation (11) need to be calibrated (or updated) due to the differences in specific physical properties between the external materials from the literature and those tested in this study. Therefore, to further improve the prediction accuracy, the model coefficients related to the G/S parameter were set to be regressed, thus yielding Equation (12).
where M r is resilient modulus; m 1 , m 2 , m 3 , m 4 , m 5 , and m 6 are model coefficients to be regressed; G/S is the gradation control parameter; σ d is deviatoric stress; θ is the first principal stress invariant; and Pa is the standard atmospheric pressure. The comparison results of laboratory-measured M r values and those predicted by Equation (12) for the above-mentioned external M r datasets are plotted in Figure 19. It shows that the prediction accuracy for the three different types of unbound aggregate materials is improved considerably, with Adj. R 2 values all greater than 0.90. Additionally, the prediction accuracy for limestone aggregate materials or those with plastic fines is relatively the best, which is consistent with the aforementioned results obtained by Equation (11). The comparison results of laboratory-measured Mr values and those predicted by Equation (12) for the above-mentioned external Mr datasets are plotted in Figure 19. It shows that the prediction accuracy for the three different types of unbound aggregate materials is improved considerably, with Adj. R 2 values all greater than 0.90. Additionally, the prediction accuracy for limestone aggregate materials or those with plastic fines is relatively the best, which is consistent with the aforementioned results obtained by Equation (11).

Summary and Conclusions
In this study, five different gradations were designed to study the gradation effect on resilient modulus behavior of TRW-derived UPAB materials. Laboratory repeated load triaxial tests were conducted under different combinations of confining pressure and cyclic deviator stress as controlled by the levels of the shear stress ratio (SSR). The Mr prediction model incorporating both gradation and stress state parameters was developed from all the 60 TRW-derived UPAB specimens tested, and the validity and applicability

Summary and Conclusions
In this study, five different gradations were designed to study the gradation effect on resilient modulus behavior of TRW-derived UPAB materials. Laboratory repeated load triaxial tests were conducted under different combinations of confining pressure and cyclic deviator stress as controlled by the levels of the shear stress ratio (SSR). The M r prediction model incorporating both gradation and stress state parameters was developed from all the 60 TRW-derived UPAB specimens tested, and the validity and applicability of the developed model were further verified by external laboratory testing database. The major findings of this study are summarized as follows.
(1) The M r values of the TRW-derived UPAB materials studied increase rapidly with increasing number of load applications. Under the same G/S value and confining pressure level, the higher the shear stress ratio (SSR) is, the greater the final M r values are, and the more significant the effect of gradation (G/S) on M r is. Under the same SSR level, the increase of confining pressure alleviates the effect of gradation (G/S) on M r of the TRW-derived UPAB materials studied. There appears to exist an optimal G/S value of around 1.6-1.8 that yields the best resilient modulus behavior of the TRW-derived UPAB materials studied. (2) The addition of standard atmospheric pressure (Pa) and the constant of 1 in the M r models has little effect on the improvement of prediction accuracy, whereas the inclusion of two additional model coefficients, k 4 and k 5 , is not effective for improving prediction accuracy of the M r models. Among the M r prediction models compared, the bivariate Model No. 8 incorporating the first principal stress invariant (θ) and deviatoric stress (σ d ) possesses the best prediction accuracy. (3) By comparing and analyzing the existing M r prediction models, a benchmark prediction model suitable for the TRW-derived UPAB materials studied was selected. The statistical correlations among the fitted model coefficients, the gradation control parameter G/S, and stress state variables were further analyzed, from which the improved M r prediction model was proposed and verified extensively to be universally applicable.