Impact of Dynamic Modulus Prediction Errors on Rutting Estimates in Sustainable Flexible Pavements

1. Introduction

One of the most common types of distress observed on flexible pavements is permanent deformation [1]. This surface distortion, known as rutting, can cause significant structural damage, reduce ride quality, and pose safety hazards due to water accumulation, reduced skid resistance, and increased risk of hydroplaning. The extent of rutting affects both the lateral and longitudinal pavement profiles, influencing drainage efficiency and vehicle stability.

Permanent deformation can result from several mechanisms, acting individually or together, including asphalt mix rutting, subgrade rutting, and densification [2,3,4]. Asphalt mix rutting occurs when depressions form along the wheel paths due to issues within the asphalt mix, such as improper compaction or deficiencies in the mix design. This distress is typically attributed to poor mix selection, inadequate design, or flaws in production. In this case, the subgrade remains stable, and the rutting is primarily confined to the asphalt layer. In contrast, subgrade rutting arises from deformations within the subgrade under repeated traffic loading. As the subgrade settles, it causes corresponding depressions on the pavement surface, compromising the structural integrity and ride quality of the road. Densification-related rutting results from insufficient compaction of individual layers during construction. Over time, as traffic loads pass over the pavement, the layers compact further, gradually leading to surface depressions and permanent deformation.

Focusing on asphalt layers, permanent deformation may consist of densification and/or flow rutting [5]. Initially, densification occurs due to compaction from traffic loads (post-compaction), while later-stage deformation results from material displacement under shear stress. Previous research [6,7] has identified three distinct rutting stages: an initial rapid densification phase, a steady-state deformation phase during the pavement’s operational lifespan, and a tertiary phase. The primary stage evolves rapidly due to post-compaction. In the next stage, rutting develops at a constant rate, mainly driven by shear stress. The final stage is characterised by accelerated plastic deformations while the material volume remains stable.

Multiple studies have explored the evolution of permanent deformation in asphalt mixtures. Common models relate permanent strain to the number of loading cycles based on a power law relationship, though they typically only capture the initial and steady-state phases without modelling tertiary deformation [8,9,10,11,12].

Given that rutting is a common distress in flexible pavements, the development of prediction models describing pavement rutting performance is a constantly evolving process. Empirical prediction models based on experience and observation lack accuracy and neglect the characteristics of the asphalt mix design [13,14,15].

In contrast, mechanistic-empirical models are more accurate in predicting rutting performance [16,17]. Moreover, machine learning (ML), including artificial neural networks (ANN) and recurrent neural networks (RNN), has been employed in the development of rutting prediction models [18,19,20,21].

In pavement design, while traditional methods set limits on subgrade strain to indirectly control permanent deformation, the Mechanistic-Empirical Pavement Design Guide (MEPDG) includes models that predict permanent deformation in asphalt pavements by considering the contributions from different pavement layers. MEPDG divides rutting contributions into four key layers: asphalt layer, granular base, subbase, and subgrade soil. Individual layer rut depths are projected over time based on traffic load repetitions [22]. The predictive models in MEPDG were initially formulated using laboratory data from repeated load tests. However, adjustments were necessary to align with field observations.

Asphalt mix properties, such as aggregate gradation, asphalt binder content, percentage of air voids, and filler, significantly influence the rutting resistance of asphalt mixtures through their impact on the stiffness modulus [23,24,25,26]. Specifically, increases in asphalt binder and air void content result in decreased rutting resistance, while an increase in stiffness modulus leads to enhanced rutting resistance [27]. Within the context of the MEPDG, asphalt mix dynamic modulus (E*) is regarded as a fundamental mechanical property of the asphalt mix. Therefore, it has a significant impact on the estimation of the rutting performance of asphalt layers during the pavement design process.

Depending on the level of analysis, MEPDG allows either the determination of E* through extensive, specialised laboratory testing or the prediction of E* using algorithms. Within the MEPDG methodology, two predictive algorithms are included: the viscosity-based NCHRP 1-37A model and the asphalt dynamic shear modulus-based model (NHCPR 1-40D). Given the importance of E* and the time-consuming procedure required to determine it in the laboratory, the need for reliable estimation through predictive algorithms is evident.

Predictive algorithms estimate E* based on mix properties, temperature, and load frequency. Numerous studies have evaluated different models, focusing on their applicability and accuracy [28,29,30,31,32,33,34,35,36]. Based on these evaluations, recalibration efforts have been undertaken, and alternative algorithms have been developed to better reflect local mixes and conditions. Although significant progress has been made in quantifying and refining the accuracy of E* predictive algorithms, limited attention has been given to quantifying and mitigating the impact of their inherent biases on pavement performance predictions, highlighting a crucial area for future research.

Previous studies have demonstrated that inaccuracies in E* prediction can lead to significant deviations in pavement performance estimation. A study by Hou et al. [37] concluded that the NCHRP 1-40D model provided better E* values than the NCHRP 1-37A model. Overall, these models were found to be useful for estimating pavement performance associated with the CAM model, particularly when actual E* tests or other performance tests are not feasible or practical in real-world applications. However, according to Georgouli et al. [38], the significant error in the E* values predicted by the NCHRP 1-37A and Hirsch models is partially mitigated in the final output, which is the fatigue cracking prediction. Still, the uncalibrated predictive algorithms led to an overestimation of fatigue cracking by 8–24%, whereas the implementation of a locally calibrated algorithm showed no significant difference.

According to Al-Tawalbeh et al. [39], analysis based on the uncalibrated Hirsch model indicated that predicted fatigue damage could vary by more than 50%, with an average deviation of 17.33%, and rutting predictions could change by up to 14%. Similarly, although the NCHRP 1-37A model demonstrated good prediction accuracy for Indiana mixtures, its use led to increased pavement distresses compared to designs based on measured E* values [40]. In contrast, limited performance deviations were observed when the NCHRP 1-40D model was applied to Egyptian mixes, reflecting the variability of model effectiveness across regions [41].

The above-mentioned studies on the effects of uncalibrated E* prediction models on pavement behaviour, especially regarding fatigue cracking and rutting, yield contradictory results, underlining the need for further investigation. Pavement sustainability best practices increasingly emphasize the importance of performance-based design supported by accurate pavement performance prediction. Reliable mechanistic–empirical modelling allows engineers to optimize pavement layer thicknesses and material selection while ensuring adequate resistance to key distresses such as rutting. This reduces premature failures and avoids unnecessary overdesign, ultimately minimizing life-cycle costs, material consumption, energy use, and environmental impacts associated with frequent rehabilitation activities. Against this background, this study examines the errors in E* estimation algorithms and their effects on rutting in asphalt pavements within the MEPDG. The tertiary rutting phase, typically associated with severe material failure, is excluded, as the adopted MEPDG-based rutting models primarily capture the initial and steady-state deformation stages. Laboratory-measured E* values (E*lab) are compared with those predicted by the NCHRP 1-37A model (E*1-37A) and a locally calibrated model developed by Georgouli et al. hereinafter referred to as the GR model (E*GR) [32]. Both measured and predicted E* values are then used as inputs in the dynamic analysis program 3D-Move [42] to simulate rutting in asphalt pavements over an analysis period. A comprehensive comparative analysis of the different rutting results is conducted.

2. Experimental Study

2.1. Asphalt Mix Properties

Multiple flexible pavement structures were considered and further analyzed (Figure 1). The pavement cross-sections differ in the thickness of the asphalt (HMA) layer, which ranges from 14 to 25 cm. The thickness of the base layer, as well as the elastic modulus of the base and subgrade layers, remain constant. The aim is to investigate whether HMA thickness, together with the estimated E* error, affects the calculated rutting.

The asphalt mix had a nominal maximum aggregate size of 19 mm and was produced using a 50/70 PEN penetration grade binder. Key mix properties, including aggregate gradation, specific gravity, asphalt content (Pb), air void percentage (Va), and effective binder content (Vbeff), were documented (Table 1).

2.2. Asphalt Mix Dynamic Modulus Testing and Estimation

Four specimens were fabricated, using a gyratory compactor, at 5% air void target. The initial dimensions of the specimens were 150mm diameter and 170mm height. The specimens were cut at the top and bottom edge and cored to meet the specified in the AASHTO T342-11 dimensions for further testing of the E*. The final dimensions of the specimens were 102±2mm average diameter and 150±2.5mm average height (Figure 2). All specimens were tested per AASHTO T342-11 at multiple temperatures (4, 15, 20, 25, 37, 54oC) and loading frequencies (25, 10, 5, 1, 0.5, 0.1 Hz). It is worthwhile to mention that although the AASHTO standard covers a wider temperature range (–10 °C to 54 °C), testing at temperatures below 4 °C was not feasible due to limitations in maintaining the required specimen temperature.

Additionally, two predictive models were used to estimate the dynamic modulus values at the aforementioned temperatures and frequencies: the viscosity-based 1-37A model [43] and GR [32] model. The viscosity-based 1-37A model was selected due to two reasons. First, it is included in input level 2 and 3 in the MEPDG software and second the global calibration factors of all asphalt concrete predictive performance models (including rut depth) have been determined using this model. The GR model was developed based on the 1-37A model following calibration to local conditions and asphalt mixes. The NCHRP 1-40D model which is also included in the input level 2 and 3 of the MEPDG software, was not utilized since according to relevant research results its predictive accuracy is significantly low for local mixes [32,38].

The general expression of both 1-37A and GR model is given in equation 1.

$$\log {\mathsf{{\rm E}}}^{\mathrm{*}}=b_1+b_2{\rho }_{200}-b_3{\left({\rho }_{200}\right)}^2-b_4{\rho }_4-b_5{V}_a-b_6\left({\displaystyle \frac{V_{beff}}{V_{beff}+V_a}}\right)+{\displaystyle \frac{b_7-b_8{\rho }_4+b_9{\rho }_{38}-b_{10}{\left({\rho }_{38}\right)}^2+b_{11}{\rho }_{34}}{1+e^{(-b_{12}-b_{13}\log \left(f\right)-b_{14}\log \left(\eta \right))}}}$$

(1)

where Ε*: dynamic modulus of mixture (psi), η: viscosity of binder (106poise), f: loading frequency (Hz), Va: air voids (% by volume), Vbeff: effective binder (% by volume), ρ34: cumulative percentage retained on 3/4 inch (or 19 mm) sieve, ρ38: cumulative percentage retained on 3/8 inch (or 9.5 mm) sieve, ρ4: cumulative percentage retained on No. 4 (or 4,75mm) sieve, ρ200: percentage passing No. 200 (or 0,075mm) sieve.

The coefficient values are shown in Table 2. Coefficient values for the 1-37A model were obtained from Andrei et al. [43], while those for GR model were derived from local calibration described in [32]. More details can be found in [38].

2.3. Rutting Model Application

MEPDG’s rutting model in the HMA layer calculates permanent deformation based on vertical plastic strain and axle load repetitions. Laboratory test results are adjusted for field conditions using calibration factors resulting in Equation 2. Parameters of the equation are explained in Table 3.

$${\Delta }_{p\left(HMA\right)}={\epsilon }_{p\left(hMA\right)}\times h_{HMA}={\beta }_{1r}\times k_z\times {\epsilon }_{r\left(HMA\right)}\times {10}^{k1r}\times n^{k2r\times \beta 2r}\times T^{k3r\times \beta 3}$$

(2)

The study performed a dynamic analysis at two vehicle speeds, 60km/h and 80km/h, using 3D-Move analysis program [42,44]. All the simulations were performed by using the typical 18-kip (8.16 ton) equivalent single-axle load (ESAL) approach. The loading configuration of the half axle is shown in Figure 3.

The wheel axle load is equal to 20kN and the pressure distributed over the circular area with a radius of 0.105m, is 0.577MPa. For simulation and comparison purposes, it is assumed that all pavement cross sections accommodate 25 million cumulative ESAL (18-kip equivalent single-axle load) repetitions over a 20-year period.

3. Data Analysis and Results

A comparative analysis assessed the predictive performance of E* models by comparing measured (E*lab) and estimated values (E*1-37A and E*GR.) (Figure 4).

The accuracy and performance of 1-37A and GR models, is described thoroughly in [38]. In brief, Root Mean Square Error (RMSE) (Table 4) calculations revealed that the GR model had the lowest deviation from laboratory-measured E*, compared to 1-37A model.

Figure 5, Figure 6 and Figure 7 illustrate the vertical strain profiles induced by a moving load traveling at 60 km/h, calculated at various depths within the asphalt layer (expressed in meters). The analysis considers dynamic modulus (E*) values obtained from laboratory testing (E*lab) and those predicted using the 1-37A (E*1-37A) and GR (E*GR) models, for total HMA layer thicknesses (h1) of 14, 20, and 25 cm. Figure 8, Figure 9 and Figure 10 present the corresponding results for a higher load speed of 80 km/h.

The developed vertical stains in every case have a similar waveform. After the main peak (~0.17s–0.18s and ~0.13.-0.14s for moving load at 60 and 80km/h, respectively), the strain response quickly dissipates, indicating that the material returns to its original state after loading. This behaviour is typical of asphalt materials, which exhibit viscoelastic recovery. At small depths near the surface vertical stress are tensile while with increasing depth until the bottom of the asphalt layer is compressive. As the depth increases, the peak strain magnitude also increases, suggesting that the material experiences greater deformation at larger depths. Furthermore, as the thickness of the HMA layer increases, the resulting strains are reduced in magnitude. A similar trend is observed with increasing load speed, higher vehicle speeds lead to lower strain levels within the asphalt layer.

Figure 11, Figure 12 and Figure 13 present a comparative analysis of the vertical strains developed at the bottom of the asphalt layer for varying HMA thickness of h1 = 0.14 m, 0.20 m, and 0.25 m, respectively and a moving load speed at 60km/h. The strains are calculated based on three different methods for estimating the dynamic modulus of the asphalt mix: laboratory testing (E*lab), the 1-37A predictive model (E*1-37A), and the GR predictive algorithm (E*GR).

According to the above figures, vertical strains tend to be overvalued when considering the E* values estimated from the 1-37A model. This is expected since E* is underestimated through the 1-37A model. However, vertical strain values calculated based on the E*GR values are very close to the ones that are obtained using the laboratory measured E* values. This is attributed to the fact that this predictive algorithm produces E* values close to the ones measured in the lab. The corresponding results for the 80 km/h load speed are omitted for brevity, as they exhibit similar trends to those discussed above.

To measure the deviations, the percentage error was computed using the vertical strain values derived from the E*lab values as a reference. Percentage error of calculated strains is presented in Figure 14, Figure 15 and Figure 16 in box-and-whisker plot. The box represents the interquartile range (IQR), which is the range between Q1 (25th percentile) and Q3 (75th percentile). A longer box suggests greater variability in the middle 50% of the data. The line inside the box indicates the median (Q2 or 50th percentile) of the dataset. The mean value is depicted with a marker. The lines extending from the box to the minimum and maximum values within 1.5 times the IQR are the whiskers. Values beyond this range are considered outliers. Any individual points outside the whiskers (Outliers), indicate potential extreme values. For clarity reasons the extreme values are not shown in the Figures below.

When considering the E*1-37A values, the percentage error in the calculated strains at the bottom of the HMA layer exhibits a wide range. For a moving load traveling at 60 km/h, the mean error is approximately −21%, −9%, and −15% for HMA layer thicknesses of 14 cm, 20 cm, and 25 cm, respectively. At a higher load speed of 80 km/h, the corresponding mean errors are approximately −18%, −23%, and 63%. Notably, for both speeds, the HMA layer with a thickness of 25 cm exhibits not only a higher mean error but also greater variability in the results, indicating increased sensitivity of the model predictions at larger pavement thicknesses.

In contrast, when the E*GR values are used, the range of percentage errors becomes considerably narrower. For a load speed of 60 km/h, the mean errors are approximately 1.9%, 3%, and 0.8% for HMA layer thicknesses of 14 cm, 20 cm, and 25 cm, respectively. Unlike the trend observed with the 1-37A model, the 25 cm HMA case exhibits both a smaller error range and a lower mean error. At the higher load speed of 80 km/h, the corresponding mean errors are approximately 10%, 10%, and −14%, with the overall range remaining nearly constant across different asphalt layer thicknesses.

It is evident that the error in the E* estimation models have a great effect on the calculated induced vertical strains at the bottom of the asphalt layer. The smaller the error the smaller the differences in the stains analysis.

Figure 17 and Figure 18 present the total permanent deformation of the asphalt layer at the end of the analysis period, calculated using E*lab, E*1-37A, and E*GR values for load speeds of 60 km/h and 80 km/h, respectively. The percentage difference of the permanent deformation is given in Table 5, with respect to the permanent deformation calculated when E*lab values are considered.

The results presented in the table indicate that the predicted permanent deformation is consistently higher when the E*1-37A values are used, with mean percentage errors ranging from approximately 24% to 27% at a load speed of 60 km/h and from 36% to 37% at 80 km/h. The differences are most pronounced at 20 cm. In contrast, the E*GR values yield substantially lower errors, within 3–4% at 60 km/h and around 9–10% at 80 km/h. These findings suggest that the GR model provides a much closer agreement with the laboratory-based E* values, demonstrating better predictive accuracy and reduced sensitivity to both load speed and HMA layer thickness. Moreover, the load speed is shown to influence the prediction errors as higher speeds correspond to larger discrepancies. This can be attributed to the viscoelastic nature of asphalt materials, where increased loading speed (or reduced loading time) leads to higher apparent stiffness and lower deformation.

The convergence of Δp(HMA) values when analysis is performed with the laboratory measured E* and those estimated with the GR model, highlights the necessity for local calibration of the E* prediction models to local conditions and mixtures.

4. Discussion

This study evaluated how the accuracy of dynamic modulus (E*) predictive algorithms affects the calculation of developed vertical strains and permanent deformation in asphalt layers of varying thickness. Using three approaches – laboratory-measured values, the 1-37A model, and the GR model – strain and permanent deformation behavior were assessed under moving load across three HMA thickness levels: 14 cm, 20 cm, and 25 cm. The main findings are summarized as follows:

• An increase in HMA layer thickness reduces vertical strains in the asphalt layer and decreases permanent deformation. This is consistent with the expected structural benefits of thicker asphalt layers, such as improved resistance to rutting.
• An increase in the loading speed (or decrease in loading time) leads to a reduction in both the vertical strains within the asphalt layer and the resulting permanent deformation. This is consistent with the viscoelastic behavior of asphalt materials, where higher loading rates cause the material to respond in a stiffer manner, reducing time-dependent deformation under moving loads.
• A detailed comparative analysis revealed significant differences in the prediction accuracy of the E* models. The GR model consistently produced dynamic modulus values that closely matched laboratory measurements, whereas the 1-37A model systematically underestimated E*.
• The agreement of the GR E* model with laboratory-determined E* values resulted in lower strain errors. The percentage error in vertical strains using E*GR values remained below 5% in all cases, with a notably narrow distribution of values, even for thicker pavements. In contrast, the 1-37A model systematically underestimated E*, leading to overestimation of vertical strains. The mean strain errors of the 1-37A model reached –21% for the 14 cm thick HMA layer, while the range of variation was considerably larger for the 25 cm thick HMA layer.
• The differences in E* and strain values are also evident in the analysis of permanent deformation. The 1-37A model overestimated Δp_HMA by more than 24% at a load speed of 60km/h and approximately 35% at 80 km/h. In contrast, predictions using the GR model exhibited only minor deviations from the laboratory results, with errors below 5% for 60 km/h and below 10% for 80 km/h.
• The thickness of the asphalt layer appears to have a greater influence on strain errors and a lesser effect on permanent deformation. It is important to note that no consistent trend was observed in strain error or rutting percentage difference with increasing pavement thickness. For example, the 1-37A model produced the highest strain error at 14 cm thickness (–21%) but also showed substantial errors at 25 cm (–15%), with a lower error at 20 cm (–9%). Similarly, the corresponding rutting percentage differences remained consistently high across all thicknesses, without a clear correlation to strain error magnitudes.
• The GR model, although more accurate overall, also exhibited fluctuations in both strain and rutting predictions without a definitive trend related to thickness. Increasing the asphalt layer thickness from 14 to 20 cm led to an increase in both strain errors and calculated permanent deformation. A further increase in asphalt layer thickness from 20 to 25 cm resulted in a decrease in both strain errors and the extent of permanent deformation.

5. Conclusions

In conclusion, the accuracy of the E* input directly affects the critical performance indicators of the pavement. Errors in predicting E* are particularly significant for pavement design regarding the rutting criterion, especially when the asphalt layer thickness is less than 20 cm. Local calibration of E* prediction models for mechanistic-empirical pavement design is therefore strongly recommended, and the choice of model is a critical factor. In the present comparative study, the 1-37A model appears to produce conservative results, while the calibrated GR model seems to optimize pavement design solutions to some extent.

Overall, it is important to have reliable predictions of pavement performance in order to achieve sustainable pavement management, since this will allow for optimal design choices, efficient resource use, and effective preventive maintenance. If pavement deterioration can be accurately predicted, then overdesign can be prevented, excessive use of asphalt resources and aggregates can be minimized, and the environmental impacts of repeated construction, such as energy use, greenhouse gas emissions, and traffic congestion, can be reduced. On the other hand, inaccurate prediction models can cause early deterioration or overly conservative designs, which will both result in increased life-cycle costs and environmental impacts. For this reason, therefore, improving the accuracy of mechanistic-empirical pavement performance predictions, especially those that are sensitive to important material properties such as dynamic modulus (E*), is an important step towards achieving smarter and more sustainable pavement infrastructure.

Future research could examine a wider range of asphalt mixes, temperatures, and traffic loading scenarios to further validate the applicability of this study’s results. Additionally, incorporating real-world performance data could enhance the robustness of predictive models and their practical integration into mechanistic-empirical pavement design systems.