Comparison of Witczak NCHRP 1-40D & Hirsh dynamic modulus models based on different binder characterization methods: a case study

The Pavement ME Design method considers the hot mix asphalt (HMA) dynamic modulus (E*) as the main mechanistic property that affects pavement performance. For the HMA, E* can be determined directly by laboratory testing (level 1) or it can be estimated using predictive equations (levels 2 and 3). Pavement-ME Design introduced the NCHRP1-40D model as the latest model for predicting E* when levels 2 or 3 HMA inputs are used. This study focused on utilizing laboratory measured E* data to compare NCHRP1-40D model with Hirsh model. This comparison included the evaluation of the binder characterization level as per Pavement ME Design and its influence on the performance of these models. E*tests were conducted in the laboratory on 25 local mixes representing different road construction projects in the kingdom of Saudi Arabia. The main tests for the mix binders were dynamic Shear Rheometer (DSR) and Brookfield Rotational Viscometer (RV). Results showed that both models with level 3 binder data produced very similar accuracy. The highest accuracy and lowest bias for both models occurred with level 3 binder data. Finally, the accuracy of prediction and level of bias for both models were found to be a function of the binder input level.


Introduction
Dynamic modulus (E*) characterizes in part the viscoelastic nature of hot mix asphalt (HMA). E* is a measure of the resistance of the HMA to deformation under compressive sinusoidal loading [1]. Several mix parameters were explored, in the NCHRP 9-19 project, in order to find which parameter strongly affects the mix performance [2].
Based on the reported results, a strong correlation was found between E* and pavement field performance such as rutting and cracking over a diverse range of traffic and climatic conditions [2].
The Mechanic-Empirical Pavement Design Guide (MEPDG), now the latest version is called Pavement ME-Design, offers three hierarchical input levels for material and traffic characterization. For the three HMA material input levels, E* is the major material property input to calculate the stress-strain relationships of the HMA layer(s) at the critical locations considering the influence of truck load, truck speed and the environmental conditions. In addition, the stiffness of the HMA layer also affects the state of stress of underlying layers.
E* values for level 1 are measured in the laboratory using the Asphalt Materials Performance Tester (AMPT) or a Universal Testing Machine. However, the E* laboratory test is time-consuming, requires costly and advanced equipment, and skilled technicians. Consequently, over years, researchers have tried to establish E* predictive models using a number of laboratory E* databases. The most well-known models are the ones developed by Professor Witczak and his students over the last few decades [3][4][5]. Hirsh model is another well-known model for E* prediction [6]. In MEPDG or Pavement ME Design, level 2 and 3 inputs for the HMA layer(s), two different E* predictive models (Witczak's models) were adopted in the Pavement ME software. The first model is called "NCHRP 1-37A" which follows a sigmoid function and predicts E* of the mix as a function of mix aggregate gradation, volumetric properties, binder viscosity and loading frequency [3]. The second model which is the NCHRP1-40D is a modified version of the latest model developed by Bari and Witczak [4]. Following is a comprehensive overview of both Witczak NCHRP1-40D and Hirsh models. Bari and Witczak [7] presented a model in 2006 to predict the HMA dynamic modulus. The model is considered as a modified version of the well-known Witczak-Andrei predictive E* model (NCHRP 1-37A model). In this model the asphalt binder is expressed in terms of Gb* DQG į QRW WKH FRQYHQWLRQDO YLVFRVLW\ Ș DV *b* can better characterize binder stiffness as a function of changes in temperature and loading frequency. Gb* is also a part of the Superpave system for the binder characterization. In 2007, this model was modified under the NCHRP 1-40D project and incorporated since version 0.95 of MEPDG under the name "NCHRP 1-40D Gb*-based E* model" [4]. The model was developed based on a large database containing 7,400 E* data points from 346 mixtures [4]. The E* database for Witczak model contains non-aged, short-term oven aged for 4 hours at 135 o C, plant aged mixes, field aged cores as well as rubber modified mixes. The mathematical structure of the model is a sigmoidal function as presented in Eq. (1)
The coefficient of determination (R 2 ) of this model is 0.91 in logarithmic scale. It should be noted that Gb* and G values contained in the database used for the model development are mostly estimated from regression models based on viscosity rather than measured in the laboratory. Thus, the error in the model's predictions may be caused by the predicted Gb* and G values. Furthermore, several studies have shown that this model significantly over predicts the dynamic modulus, particularly at the highest and lowest temperatures [8][9][10]. It was also reported that this model generally tends to overvalue the influence of temperature and underestimate the influence of other mixture characteristics such as the mix gradation [10].

Hirsh model
Christensen et al. [6] modified the original Hirsch model to predict E* of HMA. The alternate formulation of the Hirsch model "from here on will be called Hirsh Model" is shown in Eqs.
where, |E*|= dynamic modulus of the mixture, in psi; |Gb*| =|G*| = complex shear modulus of the binder, in psi; VMA= % voids in the mineral aggregates; VFA= % voids filled with asphalt; Pc= contact volume estimated from the following equation: Hirsh model was established using a much lower database containing only 206 E* measurements representing 18 HMA mixtures with only eight different binders. The R 2 of the model is 0.98. As reported by the researchers, in the above equations Gb* can be determined in the laboratory using Dynamic Shear Rheometer (DSR) or a similar device [6]. It can also be determined from mathematical models. Additionally, similar to Witczak model, it should be determined at the same temperature and loading frequency of the mixture modulus, E*, and in consistent units [6]. This model is simpler in form and number of inputs compared to Witczak model. However, it was reported to lose its prediction accuracy when used to predict E* values for Witczak database [4,7]. It was also reported by several researchers that Hirsh model, similar to Witczak model, yields significantly biased E* predictions at the extreme low and high temperatures [4,7,[11][12][13].

Comparison studies of Witczak and Hirsh models
Several studies focused on using laboratory measured E* data to compare and evaluate the accuracy of the two well-known Witczak models which are the NCHRP 1-37A and NCHRP 1-40D [14][15][16]. The studies conducted by Al-Khateeb et al. [14] and Yousefdoost et al. [15] showed that the NCHRP 1-37A model produced more accurate E* predictions compared to the NCHRP 1-40D model, particularly at the higher temperature range [. Whereas a recent  [19] concluded that Hirsh and NCHRP 1-37A models generated more accurate E* values compared to the NCHRP 1-40D model based on 1032 E* measurements from 25 mixes from the state of Idaho.
The prediction accuracy of the NCHRP 1-40D and Hirsh models was also investigated by Sakhaeifar et al at based on 3180 E* measurements. The results showed that the NCHRP 1-40D and Hirsh models produced undesirable bias at low and high E* values [20]. Recently, Georgouli et al. [21] tested 15 asphalt base and wearing course mixtures (1350 E* measurements) from Greece to evaluate the accuracy of the NCHRP 1-37A, NCHRP 1-40D, and Hirsch models. For asphalt base mixtures the results showed that the NCHRP 1-37A model produced accurate results. On the other hand, Hirsh model under-predicted E* values while Bari and Witczak model over predicted E* values. The researchers also reported that the large bias in the estimated E* values of Greece mixes using Bari and Witczak or the NCHRP 1-40D E* models is due to the associated equations used for the estimation of G* DQG į ZKLFK LQ PRVW FDVHV GLG QRW \LHOG DFFHSWDEOH YDOXHV >@ Few studies focused on the investigation of the effect of binder characterization input level, as per Pavement ME Design or MEPDG, on the prediction accuracy of the E* models [5,9,13]. Other studies evaluated the effect of the MEPDG binder input level on both Witczak E* predictive models based on 27 Idaho asphalt mixes [9,13]. MEPDG binder characterization input levels 1 and 2 were obtained from laboratory testing while level 3 binder input was based on the default MEPDG values. The researchers reported that the NCHRP 1-37A model with level 3 binder input produced the most accurate and least biased E* predictions compared to the other levels of binder input and the NCHRP 1-40D model at all three binder input levels [10]. Similar results on mixes from the Kingdom of Saudi Arabia (KSA) were also reported by Khattab et al. [5]. Thus, based on the multitude of the literature research studies covering the presented E* models, one may surmise that the prediction accuracy of these models is not consistent and vary depending of the properties of HMA and binders.

Objectives and scope of work
The main objective of this study is to compare the accuracy of Witczak NCHRP1-40D model with Hirsh model for laboratory measured E* values of 25 local KSA Superpave HMA mixtures considering the influence of the binder input level as per Pavement ME Design (level 1 versus level 2/3). Dynamic Shear Rheometer (DSR) and Brookfield Rotational Viscosity (RV) tests were conducted on the mixtures' binders using Rolling Thin Film Oven (RTFO) aged specimens. E* testing was conducted on 25 local mixes at six loading frequencies and four different temperatures. Finally, measured E* values were compared with predicted values using both The NCHRP1-40D and Hirsh models at the different binder input levels of the Pavement ME Design.

Materials
Based on the climatic regions in KSA, there are three Superpave Performance Grades (PG) that are typically used. These grades are PG 64-10, PG 70-10, and PG 76-10 which were used in this study. Eight types of modifiers (shown in Table 1) were used to achieve the required PG grades of the modified binders (PG70-10 and PG76-10). A total of 25 KSA mixes were investigated. The investigated field mixes cover the different climatic regions of KSA (i.e., central, northern, southern and, eastern regions of KSA) [5,22]. The mixtures contain14 laboratory mixes and 11 plant-produced mixes that were collected from ongoing road construction projects in KSA. Each laboratory mix had three different gradations (coarse, medium, and fine) and three percentages of air voids (2%, 4%, and 8%).

Binder testing
Laboratory tests were performed on samples from the asphalt binders contained in the investigated mixtures to determine the binder properties. Brookfield Rotational Viscosity (RV) testing, according to AASHTO TP48-97, was conducted on short term aged samples from the investigated binders using the Rolling Thin Film Oven Test (RTFO) [23]. The RV tests were conducted at three temperatures of 135, 150, and 165 °C. The Dynamic Shear Rheometer (DSR) test was also conducted on RTFO aged samples in order to determine the complex shear modulus (G*) and phase angle (į). Controls DSR model 81-PV6002 was used for the DSR testing in accordance with AASHTO T 315-06 [24]. The test was conducted at three different temperatures according to the binder PG grade and one loading frequency of 10 radians per second (1.59 Hz) in accordance with the Superpave specifications.

HMA dynamic modulus sample preparation and testing
A Gyratory Compactor was used to prepare E* samples of 150 mm in diameter and 170 mm in height. The samples were sawed from the top and bottom to obtain a final height of 150 mm and then cored in the middle to achieve a final diameter of 100 mm.
The E* tests were conducted according to AASHTO T342-11 at four different temperatures of -10, 4.4, 21.1, and 54.4 ºC. At each temperature, the test was conducted at frequencies of 25, 10, 5, 1, 0.5, and 0.1 Hz [25]. Asphalt Mixture Performance Tester (AMPT) was used for all E* tests. The majority of the E* tests were performed on two replicate samples per mix.

E* Predictions
For both E* predictive models, Gb DQG į IRU WKH ELQGHU QHHG WR EH HVWLPDWHG DW WKH VDPH WHVW temperature and frequency of E*. However, based on the Superpave requirements Gb DQG į are determined at different temperatures and only one loading frequency (1.59 Hz, which is equivalent to 10 rad/s). In order to overcome this problem, the following methodology, which is used by the Pavement ME Design, was followed to estimate Gb DQG į at the loading frequency and temperature of the measured E* values. Eq. (4) was first used to compute binder viscosity at different temperatures based on the laboratory measured Gb DQG į values at the 10 rad/sec frequency [26,27].
Finally, Eqs (6 through 11) were used to find the binder Gb DQG į IRU ERWK investigated models at the temperature and frequency at which E* needs to be determined [8]:  (11) where: fc = loading frequency in dynamic compression loading as used in the E* testing, in Hz; fs= loading frequency in dynamic shear loading mode as used in the Gb* testing, in Hz; A'= adjusted "A" (adjusted for loading frequency); VTS'= adjusted "VTS" (adjusted for loading IUHTXHQF\ Ș ELQGHU YLVFRVLW\ DV D IXQFWLRQ RI ERWK ORDGLQJ IUHTXHQF\ Is) and temperature (TR), in cP; TR = temperature, in Rankine; |Gb*|= Gb* = complex binder shear modulus, in Pa.

Pavement ME-design binder input levels
For input level 1a (Pavement ME Design level 1 conventional binder data), the RV test results on the RTFO aged binders were fit with Eq. (5) to find the A-VTS values. For input level 1b (Pavement ME Design level 1 Superpave performance grade (PG) binder data), the DSR results on the RTFO aged binders were used to calculate the viscosity at each E* tested temperature using Eq. (4). The A and VTS values for each binder were calculated using Eq. (5). Finally, Equations (6 through 11) were used to estimate Gb DQG į DW WKH VDPH frequency and test temperature of E*. Lastly, for input level 3 (Pavement ME Design level 3 default A-VTS parameters), the Pavement ME Design default A-VTS values stored in the software were selected based on the Superpave binder PG grade. Table 1 summarizes the A and VTS parameters estimated for levels 1a and 1b. The coefficient of determination (R 2 ) as well as the default A and VTS values (level 3) recommended by the Pavement ME Design for the investigated binders are also shown in the table. The high R 2 values indicate excellent fit.
The bias in prediction was assessed by determining the closeness of the slope and intercept of the regression line of measured and predicted E* values to one and zero, respectively [29]. The effect of the binder input level on the predicted E* values using the

ASCMCES-17
7003 NCHRP1-40D model is shown in Figure 1 (a through c). The figure shows relatively very good prediction accuracy as indicated by the goodness of fit statistics summarized in Table  2. The prediction accuracy (lowest scatter) was the highest for the predictions based on level 3 binder characterization (Pavement ME Design default A-VTS values based on the PG grade of the binder) with R 2 of 0.84 and Se/Sy of 0.40. However, all three cases showed variable degrees of bias in the predictions as indicated by the slope and intercept of the regression line shown in the figures. Again, the lowest scatter occurred when level 3 binder values was used (slope of 0.9402 and intercept of 0.365). The unconstrained regression lines shown in Figure 1 (a to c) also show that the bias in the predictions based on level 1a binder data was mostly in all E* values while for the case of level 1b and level 3 binder data it was mostly in the low and intermediate E* values.
For Hirsh model the comparison of measured and predicted E* values at the different binder characterization levels is shown in Figure 1 (d to f). The goodness of fit statistics for the E* predictions shown in Figure 1 indicates a very good accuracy (low scatter). Similar to Witczak model, the data suggests that the highest accuracy (R 2 = 0.85 and Se/Sy = 0.39) and lowest bias (slope = 0.862 and intercept = 0.4107) occurred when level 3 binder data was used. For Hirsh model the bias in the predictions was mostly at the high and low temperatures (low and high E* values) as reported by many research studies [4,7].
The figure shows that the accuracy of both models when E* is predicted based on either level 1b or level 3 binder is very close. If level 1a binder data is used, then Hirsh model is more accurate compared to Witczak model. However, the bias of Hirsh model is higher compared to Witczak model at all binder input levels especially at the high and low E* values.

Conclusions
The accuracy and level of bias of both Witczak NCHRP1-40D and Hirsh dynamic modulus prediction models with three Pavement ME Design binder characterization input levels were studied for 25 local KSA Superpave mixtures. The results of this study suggest that the performance (in terms of accuracy and bias) of the investigated models is greatly influenced by the binder characterization method. Both Witczak NCHRP 1-40D and Hirsh models yielded relatively accurate E* predictions for the KSA mixes with very similar goodness of fit statistics for the predictions based on level 1b and level 3 binder data. However, Hirsh model along with binder data based on level 1a (A-VTS values based on RV data) was more accurate compared to the NCHRP 1-40D model. Both models showed some bias in the predictions which varied with the binder input data. For Hirsh, the bias was more towards the lower and higher E* values while for Witczak the bias was more towards the lower and intermediate E* values. Generally, the NCHRP1-40 model yielded lower bias compared to Hirsh model at the three binder characterization input levels. Thus, Hirsh model needs some improvements in terms of bias. Both Witczak and Hirsh models with level 3 binder data produced very similar accuracy based on the goodness of fit statistics. Finally, the highest accuracy and lowest bias for both models occurred with level 3 binder data.