Abstract
This study presents a new mechanistic procedure for determining the critical cracking temperature of asphalt concrete (AC) using data from bending beam rheometer (BBR) test of asphalt binder and indirect tension (IDT) test of AC. This new procedure uses BBR creep data to generate the mixture relaxation modulus mastercurve by utilizing the Hirsch model, time-temperature superposition principle, and Prony series-based interconversion method. The Hirsch model parameters are calibrated by comparing creep data from BBR and IDT creep tests performed at the same temperature. Boltzmann hereditary integral and second-order heat equation are then used to calculate thermal stress from the developed relaxation modulus mastercurve. IDT strength data is transferred from test strain rate to thermal strain rate using the viscoelastic continuum damage model. Since a strain gauge is not attached for traditional laboratory IDT strength testing, this study derived an analytical equation based on the Hondros solution to compute the horizontal strain rate from the applied vertical displacement rate. Finally, the critical cracking temperature is determined by coupling the thermal stress and strength profiles. Using the procedure presented in this paper, the critical cracking temperatures of four AC mixtures were predicted from BBR and IDT data. Their actual critical cracking temperatures were measured using thermal stress restrained specimen test performed in the laboratory to validate the method. The predicted critical cracking temperatures are found to be very close to the laboratory measured values. The developed procedure has substantial practical and technical importance in predicting the critical cracking temperature of AC because it utilizes widely available BBR and IDT tests.
Low-temperature cracking is one of the major distresses that can significantly reduce the service life of an asphalt concrete (AC) pavement. This type of crack mostly occurs in colder regions where the temperature drops significantly below the freezing point. It can be initiated by a single extreme low-temperature event or by thermal fatigue cycles (1, 2). As the temperature decreases, the pavement layer tries to shrink its volume because of the thermal contraction of the material. Since the thermal strain is unreleased, tensile stress develops within the pavement layer. When the thermal stress exceeds the tensile strength of the AC layer, micro-cracks are initiated. The temperature at which the thermal stress exceeds the tensile strength of AC is known as the critical cracking (or failure) temperature ( 3 ). After formation, these micro-cracks penetrate the full depth and expand to the entire width of the AC layer depending on the temperature cycles and the traffic loadings. The low-temperature cracks are observed mostly in the transverse direction of the pavement and are alternatively known as transverse cracks. They can be a source of further structural deterioration and can lead to premature failure of a pavement. Since weather conditions vary significantly from location to location ( 4 ), it is necessary to select a proper asphalt mixture that can withstand against the cracking at the lowest temperature during its lifetime. In other words, the critical cracking temperature of an AC mix should be known before using it to ensure adequate service life.
The thermal stress restrained specimen test (TSRST) can be very helpful to determine the critical cracking temperature of an AC mix ( 5 ). In this test, a prismatic AC specimen is restrained at both ends from movement and a continuous temperature drop is applied until the specimen fails. It is a temperature-controlled test that determines the thermal stress profile and failure temperature of an AC specimen. Though the test process is similar to the field low-temperature cracking event, it is still not widely accepted from some practical points of view. The main problem is that the TSRST is an expensive test because it needs extremely low temperatures to complete a test. The ordinary environmental chambers that are available in most of the AC mix design laboratories are not suitable for this test. In addition to cost, the sample ends need to be attached to two end platens with a special kind of glue that must be capable of handling extremely cold conditions during stress measurement using a loading frame and controllers. The process of sample alignment with end platens is extremely difficult, and the repeatability of the test is diminished by improper alignment ( 6 ). Furthermore, the actual time required for completing one test is unknown. For these reasons, researchers are exploring alternative approaches to determine the critical cracking temperature of AC mix.
As an alternative to direct measurement, an analytical approach can be useful to predict the critical cracking temperature of an AC mix. In the analytical approach, thermal stress needs to be calculated first from the material’s stiffness and rate of temperature drop. Then, it can be coupled with the strength profile to determine the critical cracking temperature. However, asphalt binder/mixture is a viscoelastic material. The prediction of its thermal stress profile is not as simple as it is in some other elastic materials, because it is highly temperature susceptible material and its response to a load depends not only on the magnitude of the load but also on the time of loading. It is necessary to apply the Boltzmann hereditary integral form of the constitutive equation of a linear viscoelastic material to relate time-dependent stress and strain ( 7 ). This approach uses the relaxation modulus of an AC sample to develop its thermal stress profile for a constant rate of temperature drop. The current standard test to determine the tensile creep compliance (which can be transferred into relaxation modulus) of an AC sample is the indirect tensile (IDT) test (AASHTO T 322). Besides the creep property, this test is also used to determine the IDT strength of an AC sample. There are several reasons to choose the IDT test over the uniaxial tension test: for example, it is a simple test procedure, the failure plane is known a priori, and state of stress is similar to the state induced under wheel loads ( 8 ). According to the current approach, the IDT creep test is performed on the AC specimen at different temperatures to determine its creep performance. The creep data is then used to develop the thermal stress curve for a constant rate of temperature drop. Next, the IDT strength test is performed at different temperatures to determine the tensile strength profile. Finally, the critical cracking temperature is determined by coupling the thermal stress curve with the IDT strength curve ( 9 ). However, the stress-strain relationship of a viscoelastic material depends not only on the temperature but also on the applied stress/strain rate. The IDT strength test is done by applying a fast vertical crosshead displacement rate (12.5 mm/min) on the sample, which is higher than the thermal contraction rate. Thus, the direct comparison between the strength data from the IDT test and thermal stress computed from the creep data is fundamentally erroneous. On the other hand, performing IDT strength test at thermal contraction rate is not feasible. Moreover, it is not recommended to attach strain gauge fixtures while performing the IDT strength test ( 10 ), because failure in the low-temperature IDT test is abrupt and may destroy the strain gauge. Thus, the horizonal strain is not recorded while performing the IDT strength test while it is recorded during the IDT creep test. None of the previous studies reported any analytical approach that can compute horizontal strain rate from the crosshead displacement rate. Furthermore, thermal stress was computed for the uniaxial condition, whereas the stress-strain condition in the IDT test is biaxial ( 11 ). Additionally, the IDT test at low temperature needs a special environmental chamber which is not available to most of the asphalt production laboratories.
As testing of an AC sample at extremely low temperature is difficult, it would be very useful if traditional asphalt binder test data, such as from bending beam rheometer (BBR) and direct tension (DT) testing, could be used to determine the critical cracking temperature. Bouldin et al. ( 3 ) proposed a method that computes the thermal stresses using BBR data. The binder thermal stresses are then converted to mixture thermal stresses by multiplying with an empirical constant equal to 18. Next, the DT test is used to determine the tensile strength of the binder at different temperatures. Finally, the critical cracking temperature is predicted by comparing thermal stress and strength profiles. However, this method has two fundamental limitations. The first limitation is related to the empirical pavement constant because mixture properties depend not only on binder properties but also on mixture volumetrics. Thus, use of single pavement constant is not enough for all types of mixtures. Roy and Hesp ( 12 ) estimated that the value of pavement constant can vary from 3.4 to 16.7, which is a relatively large range. Thus, use of a single pavement constant value may give significant errors in predicting critical cracking temperature. Some researchers have reported that the Hirsch model for converting binder data into mixture data can be used to eliminate the empirical pavement constant (13, 14). The other problem of the method of Bouldin et al. ( 3 ) is that the strain rate (3% per min) applied to measure tensile strength in the DT test is very high compared with the thermal contraction rate. The repeatability of DT testing is also an issue. As such, no accurate mechanistic method is available yet to combine the thermal stress profile and strength data to predict the critical cracking temperature of an AC mix.
As thermal cracking reduces the service life of a pavement significantly, proper prediction of the critical cracking temperature is crucial for ensuring a long-lasting pavement. This study develops a mechanistic procedure to predict the critical cracking temperature of an AC mix accurately without performing TSRST. To do so, this study uses data from conventional BBR and IDT tests. It can be noted that the BBR is a widely used test on asphalt binder to capture its low-temperature creep properties and performance grade (PG). The IDT test is extensively used during the mix design stage to determine the indirect tensile strength, but it can also be applied to capture the viscoelastic properties of an AC mixture. Unlike the previous works, this study includes a mechanistic procedure to transfer the IDT strength data from test strain rate to thermal contraction rate. The developed procedure was applied to predict the critical cracking temperatures of four different AC mixtures, made out of four different types of binders. The procedure is validated by comparing the predicted values with the laboratory measured TSRST values.
Development of the Method
The method developed in this study has three primary modules: development of mixture relaxation modulus from the binder test data, calculation of thermal stress occurring with constant rate of temperature drop, and conversion of IDT test strength from test strain rate to thermal contraction rate. Details of these modules are given below.
Mixture Relaxation Modulus from BBR Test
In the BBR test, an asphalt beam sample is subjected to flexural loading over a total time of 240 s to capture the creep behavior of the binder at any test temperature. As a viscoelastic material, the creep deflection occurs in the beam over the time. The applied load in the BBR test is significantly small, and it is assumed that no plastic bending occurs in the beam during the entire test. Therefore, the classical beam theory, which assumes the linear distribution of stress as well as strain along its depth, can be applicable to the BBR test data to calculate the flexural creep stiffness of the asphalt beam. The flexural creep stiffness is the ratio obtained by dividing the maximum bending stress in the beam by the maximum bending strain. Both maximum bending stress and strain occur at the top and bottom of the beam at its midspan. They can be computed from the applied load, creep deflection at the middle of the beam, and beam geometric information using classical beam theory. The final expression of the flexural creep stiffness, S(t) is given below:
where P is the applied load; δ(t) is creep deflection at the middle of the beam; l, b, and h are the length, width, and depth of the sample; t represents the time. It should be noted that the flexural creep stiffness S(t) is the reciprocal of creep compliance D(t) and a numerical interconversion method is required to transfer it into relaxation modulus E(t). It is known that the BBR test is performed at different temperatures to capture the temperature susceptibility of the binder. The time-temperature superposition principle (TTSP) can be used to construct the creep stiffness master curve with respect to a reference temperature. For fitting the master curve, this study used the CAM model ( 15 ) as shown in Equation 2.
where Sg is the glassy modulus that can be determined by fitting the test data; β, κ, and λ are the fitting parameters that also can be found by fitting tested data; and ξ is the reduced time that can be calculated from the loading time t and shift factor aT using Equation 3.
The value of aT depends on the test temperature and the reference temperature at which the mastercurve is drawn. Currently, several different equations are available for determining aT. However, this study used the second-order polynomial expression, as shown in Equation 4, for fitting aT to make this analysis simple.
where T is the test temperature (°C); C1, C2, and C3 are the fitting parameters. After constructing the creep stiffness master curve, the Prony series-based interconversion method ( 16 ) can be used to convert the creep stiffness master curve into the relaxation modulus master curve. To do so, the creep stiffness master curve of the binder is first transferred into the creep compliance master curve D(t) by taking the reciprocal values and the compliance data is fitted to a generalized Kelvin model, as shown in Equation 5.
where Dg is the glassy compliance, Dj is the retardation strengths, and τj is the retardation times. Similarly, the unknown relaxation modulus can be represented by the generalized Maxwell model, as shown in Equation 6.
where Ee is the equilibrium modulus, Ei is the relaxation strengths, and ρi is the relaxation times. The relationship between creep compliance and relaxation modulus can be expressed by Equation 7:
Using this relationship, the unknown values of Ee and Ei can be determined by solving a system of equations. Before that, the unknown ρi values need to be determined from the graphical root finding method after converting Equation 7 into Laplace domain. Using the values of Ee and Ei, the binder relaxation modulus E(t) can be determined from Equation 6.
After determining the binder relaxation modulus, the mixture relaxation modulus can be determined by applying an appropriate composite material model. Currently, there are a few models ( 17 – 20 ) available for predicting mixture dynamic modulus from binder modulus and mixture volumetric properties. These predictive models can be divided into two categories: statistical models (e.g., Witczak models) and semi-empirical models (e.g., Hirsch or Alkhateeb models). The statistical models are developed based on statistical relationships of mixture modulus with binder modulus and mixture volumetric properties. Large amounts of data are needed to calibrate the model parameters. On the other hand, the semi-empirical models are developed based on the law of mixtures. They can predict the mixture modulus by employing a three-phase system of aggregate, binder, and air void in series, parallel, and/or combined arrangement. Several past studies evaluated the effectiveness of these predictive models and they concluded that neither of these models can provide accurate prediction without proper calibration ( 21 – 23 ). It is also recommended to calibrate the model parameters if the loading mode is different than the loading mode in the dynamic modulus test ( 24 ). This study chose the semi-empirical Hirsch model over statistical models to transfer binder relaxation modulus into mixture relaxation modulus. Moreover, Zofka et al. ( 25 ) also reported that the calibrated Hirsch model can be used successfully to interconvert mixture and binder low-temperature properties. The expression for the Hirsch model is given below:
where Emix(t) is the mixture relaxation modulus at time t in MPa; Eagg is the modulus of the aggregate which is assumed to be insensitive of time and temperature; and Eb(t) is the binder relaxation modulus at time t in MPa; VMA is the percentage of void in the mineral aggregate in the AC mixture; VFA is the percentage of void filled with asphalt in the AC mixture; and Pc is an empirical contact factor that can be determined from Equation 9.
where H1, H2, and H3 are the regression parameters, and they can be determined by comparing the binder and mixture creep test data performed at same temperature. After deriving the mixture relaxation modulus from the Hirsch model, it is fitted into Equation 2 to determine the mastercurve parameters using the previously calculated reduce time ξ.
Calculation of Thermal Stress from Mixture Relaxation Modulus
The thermal stress σ(t) in a viscoelastic material can be determined using the Boltzmann hereditary integral after considering an imaginary thermal strain, as shown in Equation 10.
where E(t,T) represents the relaxation modulus of the mixture at time t and at temperature T; ε is the thermal strain; and τ is an arbitrary time variable. In Equation 10, the relaxation modulus E(t,T) is a function of both time and temperature, and the TTSP can be used to convert it to a function of time only. On the other hand, the thermal strain ε is a function temperature. If the coefficient of thermal contraction CTC is assumed to be constant, the rate of thermal strain can be expressed as follows:
Since the temperature T is a function of time, the rate of thermal strain dε/dt will also be a function of time. After plugging Equations 2, 3, 4, and 11 into Equation 10, the final thermal stress equation can be written as follows:
where Eg, β, and κ are the coefficients of the mixture relaxation modulus mastercurve fitted into Equation 2; and C1, C2, and C3 are the shift factor coefficients.
Since the validation of the proposed method needs to be done by comparing predicted data with actual TSRST results, temperature distribution inside the specimen should be considered. Figure 1 represents a schematic of the TSRST test. In this test, a prismatic (cylindrical) AC test specimen is glued to two end platens. The glued specimen is then placed in an environmental chamber of the load frame in such way that one end is attached with fixed support and other end is attached with a load cell, as shown in Figure 1a. The temperature of the chamber is then cooled down at a constant rate. During the test, both ends of the specimen are kept in constant positions to prevent deformation of the specimen. As the temperature decreases, the thermally induced stress gradually increases, as shown in Figure 1b. Initially, the stress increment rate is slow because of the higher relaxation property of the mixture. As the temperature decreases further, stress increases rapidly until the sample fails. Since it is not possible to measure the stress directly, the thermal load over time is recorded by the load cell and stress is computed by dividing the load by the cross-sectional area. Thermocouple is used to measure the temperature outside the sample, inside a dummy sample, or both. The temperature distribution and consequently the thermal stress distribution in a cylindrical sample are not linear, as shown in Figure 1c. Thus, use of only outside or inside temperature (or their arithmetic mean) will produce error in further analysis. The effective mean temperature T(t) of the sample at time t can be determined by applying integration over the radius, as shown in Figure 1d:
where r0 represents the radius of the sample; and T(r,t) represents the temperature at any location r from the center at time t. The temperature profile T(r,t) can be determined by solving the second-order heat equation in polar coordinates by considering one-way heat transfer through the radial direction:
where αT is the thermal diffusivity of the sample. If the initial temperature condition Ti and temperature drop rate C0 are known, then T(r,t) can be determined by solving the above equation. The solution is shown in Equation 15 ( 26 ).

Schematic of thermal stress restrained specimen test (TSRST): (a) test setup and loading condition, (b) temperature and stress over time, (c) longitudinal sectional profile, and (d) cross-sectional profile.
where F0 is known as the Fourier number and F0=αT t/r0 2 ; J0 and J1 are the Bessel functions of the first kind and can be computed using Equation 16; μ represents the roots of the Bessel function J0 and can also be found elsewhere ( 27 ).
This study found that first five terms of the infinite series of Equations 15 and 16 can give reasonable accuracy, thus higher order terms can be dropped. For further computational simplification, the effective mean temperature can be fitted into a polynomial equation and combined with Equation 12 to determine the thermal stress as a function of time. It is not possible to solve the integration analytically, but it is possible by numerical approximation.
Conversion of IDT Strength into Thermal Strain Rate
The IDT strength test is used to determine the IDT strength of the AC by applying a constant crosshead displacement rate on the sample along the vertical axis. Figure 2 shows a schematic of the IDT test. Figure 2a represents the loading condition in the IDT test. Hondros ( 28 ) provided an analytical solution, considering it as a plane stress problem. Based on this solution, the response stress distributions along vertical and horizontal axes are shown in Figure 2b. If an IDT specimen with radius R and thickness d is subjected to a strip load of width a = 2R sin(α) and magnitude P, the horizontal and vertical stresses at any point y on the y-axis can be found by Equations 17a and 17b, respectively.

Schematic of indirect tensile (IDT) test: (a) loading condition and (b) stress distribution profiles.
As the loading condition in IDT test is biaxial, the horizontal and vertical strains can be determined after using Hooke’s law in the biaxial plane if the modulus E and Poisson ratio ν are known:
Because of the nature of the IDT specimen geometry, another constant can be derived to correlate the horizontal and vertical strains:
Like strain, the vertical displacement can be computed by integrating the vertical strain, as shown in Equation 20.
If the applied vertical crosshead displacement of the top point uy(0,R) is known, then the horizontal strain at any point εx(0,y) can be determined after combining Equations 18a and 20:
Again, for the plane stress condition, the stress can be written as follows:
After substituting Equation 19 into Equation 22, the horizontal stress can be computed as a function of horizonal strain.
In an IDT test, the maximum horizontal stress occurs at the center of the specimen (y = 0). Theoretically, the maximum horizontal tensile strain should be at the center of the specimen and the failure plane should be the vertical diametral plane. However, AC is composed of aggregates, asphalt mastic, and voids. Because of the variation in stiffness of these materials, the strain distribution in such a composite material system differs slightly from the theoretical analysis. Zhang et al. ( 29 ) reported that AC can be regarded as a homogenous material only on a representative volume basis. The creep test is performed by measuring deformation over a sufficient gauge length. To be consistent with the creep test, this study prefers to use the average values at middle of the specimen instead of extreme values at just the center of the specimen. The average values can be determined as follows:
where l represents half of gauge length and was assumed to be 19 mm (the standard gauge length is 38 mm for a 150 mm diameter sample). For simplicity, the “avg” term in the subscript of each expression is omitted in the following parts of this paper. The expression for the horizontal stress can be written as:
Equation 25 is derived for an isotropic elastic material, but it can be modified for a linear viscoelastic material by replacing viscoelastic counterparts. It should be reported that, like stiffness, the Poisson ratio of a viscoelastic material is also a function of time. However, this study used a constant Poisson ratio instead of a time-dependent value to make the analysis simple.
where horizontal strain rate can be determined by differentiating Equation 21 with respect to time. It should be noted that the machine compliance is ignored as the stiffness of steel is significantly greater than the stiffness of AC.
After obtaining the horizontal strain rate, it is necessary to transfer the tensile strength of the sample from the IDT test strain rate to the thermal contraction rate as well as from the test temperature to other temperatures. This study prefers to use the viscoelastic continuum damage (VECD) model, because past studies showed that it is valid for any loading mode and the TTSP can be used with it ( 29 – 32 ). The only concern for the VECD model is that the damage in the sample must be caused by the accumulation of micro-cracks and this is true for the IDT test. It should be noted that the damage in the IDT test is the result of both compression and tension effects. However, they are linked with each other and this study considered only tensile effect for the further analysis. The VECD model uses the viscoelastic corresponding principle (stress-pseudostrain), along with the work potential theory proposed by Schapery ( 33 ), to derive the relationship between pseudo-stiffness C and the internal state damage variable S of a sample. From the viscoelastic corresponding principle, the stress σ at any time t in a damaged sample can be calculated using Equation 28.
where C is the pseudo-stiffness, ER is an arbitrary modulus (assumed to be 1), and εR is the pseudo-strain calculated from the relaxation modulus E(t) and applied strain rate dε/dt for the uniaxial condition:
However, for the IDT test condition, the expression for the pseudo-strain will be changed according to Equation 26.
From the potential work theory, the damage evolution law is given below:
where WR is the pseudo-strain energy defined as the area under the stress versus pseudo-strain curve and can be computed by using Equation 32. S is the internal state damage parameter. α is the damage evolution rate that can be obtained as (1+1/m), in which m is the linear slope between stiffness and frequency in a logarithmic scale.
After substituting Equation 32 into Equation 31, the internal state damage variable S can be computed from the test data using chain-rule after considering C as a function of time, t:
where I is a parameter used for minimizing the variability between samples. As mentioned earlier, the relationship between C and S is a unique property of a material and this relationship is independent of loading type, rate, and mode. It can be expressed as follows:
where C11 and C12 are the material parameters. It should be noted that, depending on the applied load, the C at failure changes and is always greater than zero. The characteristic damage curve in Equation 34 is developed from the IDT test data, but it can be used to predict stress-strain for the uniaxial load condition by using pseudo-strain from Equation 29. Now, Equations 28, 29, 31, and 32 can be combined to give Equation 35, which correlates the internal damage of the material S to the strain history of the material. After integration of Equation 35 with respect to time, this relationship can be used to compute the value of S for any loading condition.
The found value will be used to calculate the C and so σ(t) using Equations 34 and 28, respectively. Finally, the tensile strength can be determined by differentiating Equation 28 with respect to time. It will give the rate of stress development with regard to time.
In Equation 36, the left-hand side represents the rate of stress. On the right-hand side, the first term is the product of the pseudo-stiffness rate and the pseudo-strain, and the second term is the product of the pseudo-stiffness and the pseudo-strain rate. In monotonic testing, damage always increases as load increases. Therefore, the pseudo-stiffness rate will always be negative, and the pseudo-strain will always be positive. Thus, the first term on the right-hand side is always negative. On the other hand, the pseudo-stiffness is always equal to or greater than zero, and the pseudo-strain rate is obviously greater than zero. Thus, the second term on the right-hand side of the equation is always positive. Now, the right-hand side defines a curve that starts from a positive value, and by passing from a point that corresponds to the strength value it enters a negative region. The strength can be assumed to be the maximum stress when this rate of stress will be zero. The abovementioned VECD method is derived to find the strength for any loading rate at one specific temperature. Daniel and Kim ( 31 ) showed that TTSP can be used with the VECD model to predict stress-strain in damaged linear viscoelastic material. This can be done by replacing the time variable t with the reduced time ξ from Equation 3 in the respective equations. In this way, the strength of an AC mixture at thermal strain rate can be found for different temperatures. The literature reports that IDT strength increases as the temperature decreases up to a maximum value, then it decreases as the temperature further decreases ( 34 ). The temperature at which the IDT strength reaches its maximum values is around –10°C. Differential thermal contractions between aggregate and binder at further lower temperatures cause internal damage within the sample and reduce the IDT strength. However, it is difficult to determine the internal damage and it is beyond the objective of this study. Thus, the characteristic damage curve at –10°C is assumed to be the most favorable for the sample, and it can be used to compute the AC mixture’s strength at further lower temperatures. Finally, the critical temperature can be found by locating the intersection of thermal stress and strength curves.
Summary of the Proposed Procedure
The proposed procedure has five primary steps, which also have their sub-steps, as shown in Figure 3. The outline is as follow:
Step 1: To calibrate the Hirsch model from BBR and IDT creep test data.
Step 2: To develop the relaxation modulus master curve of the AC mixture from BBR test data using TTSP, Prony series-based interconversion method, and Hirsch model.
Step 3: To predict the effective temperature of the sample using the second-order heat equation, and to calculate thermal stress from the mixture relaxation modulus using Boltzmann hereditary integral.
Step 4: To determine the horizontal strain rate from vertical crosshead displacement rate using Hondros solution, and to convert IDT strength from test strain rate to thermal contraction rate using the VECD model.
Step 5: To determine the critical cracking temperature by locating the intersection of the thermal stress and tensile strength curves.

Summary of the proposed procedure.
Practical Implementation
The mechanistic procedure developed in this study is intended to determine the critical cracking temperature of an AC sample without performing TSRST. Though this procedure consists of several numerical methods, most of them are already being used in the pavement industry. A computer program can easily be developed to perform these numerical methods. The end user just needs to input the test data and the program will give the critical cracking temperature value by executing the numerical methods. Thus, the developed procedure can be performed in any normal pavement engineering laboratory. Moreover, this procedure does not require any new test or any new equipment, rather it uses data from the widely available BBR and IDT tests. Thus, the proposed method will also be an economical solution for determining the critical cracking temperature of an AC mixture.
Materials and Tests
To apply and validate the proposed method, this study used four different AC mixtures, made out of four different types of binders—one unmodified and three SBS (styrene-butadiene-styrene) polymer modified—at three different concentrations (1%, 3%, and 5%). The BBR test was performed on the binder samples to determine their low-temperature creep behavior. The test was done at four different test temperatures (0°C, –6°C, –12°C, and –18°C). Three replicate samples were used for the BBR test. From this test, the creep deflection values were recorded at six loading times: 8, 15, 30, 60, 120, and 240 s. The BBR test was performed on long-term aged binders. The binder creep data at 0°C were compared with the mixture creep data to calibrate the Hirsch model parameters. The complete BBR test data was then used to develop the mixtures’ relaxation modulus master curves and thus the mixtures’ thermal stress curves using Boltzmann integral. For the AC mixture samples, a binder content of 4.7% was determined using the Superpave mix design procedure for the base (unmodified) binder for an SP-III (gradation) aggregate. The binder content and aggregate gradation were kept the same for the modified binders to make the study simple and consistent. Similarly, the researchers aimed to keep the targeted air void (7.0±0.3%) the same for all the tested mixture samples. The details of aggregate gradation and mix volumetric are reported in Table 1.
Aggregation Gradation and Mix Volumetric Properties
This study used cylindrical AC samples with a diameter of 150 mm and a height of 38 mm for the IDT test. At the beginning, the IDT creep test was performed on each specimen at 0°C to determine the creep compliance by following the procedure specified in AASHTO T322. After the test, the specimen was rested for 1 h for recovery. Finally, the strength test was performed (also at 0°C) by applying a constant crosshead vertical displacement rate of 12.5 mm/min (0.208 mm/s). The response load over the time and failure load of each sample was recorded. Then, Hondros solutions were applied to compute the horizonal stress in the middle of the sample. The proposed Equation 27 was used convert the vertical crosshead displacement rate into the horizontal strain rate. Later, the VECD model was used to transfer the IDT test strength from test strain rate to thermal contraction rate. After that, the critical cracking temperature of each mixture was determined from the intersection of thermal stress and strength curves. The measured thermal diffusivity and coefficient of thermal contraction of the mixture samples are 4.0 × 10–7 m2/sec and 2.0 × 10–5 mm/mm/°C, respectively. The Poisson ratio of mixtures was determined as 0.20. For validation, this study performed the TSRST test on cylindrical AC samples with a diameter of 50 mm and a height of 170 mm. Before the test, each sample was glued with two end platens by plastic steel putty and cured for more than 24 h to ensure sufficient bonding. The sample was then placed in the environmental chamber at 0°C for a sufficient amount of time to achieve thermal equilibrium. Finally, the test was performed by applying a cooling rate of 10°C/h (2.8 × 10–3°C/s) until the sample failed. From this test, the thermal load over the time, failure load, and critical cracking temperature of each sample were recorded. The thermal stress was calculated by dividing thermal load with the cross-sectional area of the specimen. For simplicity, this calculated stress is referred as the measured thermal stress in this paper. Finally, the predicted and measured values were compared to evaluate the effectiveness of the proposed method. In this study, two replicate samples were used for all mixture tests and the average values are reported in this paper.
Results and Discussion
At the beginning, the Hirsch model parameters were calibrated by comparing the creep data from the binder BBR and AC mixture IDT creep tests, performed at the same temperature (0°C). To gain better accuracy, the model was calibrated separately for each binder. First, Equation 8 was rearranged to calculate the contact factor Pc:
To make the calibration process easy, instead of the relaxation modulus, the creep stiffness was directly used in the above equation to determine Pc. It should be noted that both binder and mixture creep stiffnesses must be measured at the same loading times (and also at the same temperature) to work this procedure appropriately. Since the binder creep stiffness was measured at six loading times, the number of Pc values was also six. Next, the found Pc values and other inputs were used in Equation 9 to determine the Hirsch model parameters (H1, H2, and H3). The Microsoft Excel solver was applied in this study for this calibration process. Table 2 lists the calibration parameters found for each binder type. For all cases, good fittings were observed as the found coefficient of determination, R2 values, were close to unity.
Calibrated Hirsch Model Parameters
In the next step, TTSP was applied to construct the creep stiffness mastercurve of each binder from the BBR data performed at four different temperatures. The reference temperature of the constructed master curve was selected as 0°C. Next, the creep stiffness values were used to calculate the values of Pc using Equation 9. Then, the calculated Pc values were used to transfer the binder creep stiffness into mixture creep stiffness using Equation 8. After that, the mixture creep stiffness was converted to creep compliance by taking the reciprocal values. Finally, the mixture relaxation modulus was determined from the mixture creep compliance using the Prony series-based interconversion method. To fit Prony series to the laboratory test data, a dataset of discrete-times equally spaced in the logarithmic scale was created, where the smallest time was a few decades smaller than the lowest reduced time of the creep compliance mastercurve and the largest time was a few decades larger than the largest reduced time. Since this dataset had a wide range of discrete-times, it could fit the test data almost perfectly and the sum of the squared error was found to be very small. Next, the lowest discrete time was dropped from the discrete time dataset, and the new discrete time dataset was used to fit the test data. Then, the new sum of the squared error was calculated and compared with the previous value. If the value was stable (almost the same), then the next lowest point was dropped, and this process was continued until the sum of the squared error showed an abrupt change. The time at which the sum of the squared error remained stable was considered as the optimum lowest discrete time (minimum). In a similar way, the optimum largest discrete time (maximum) was determined. Finally, the optimum number of Prony terms was determined by counting the discrete time between the optimum lowest and optimum largest discrete times. The optimum number of Prony terms was found to be five. Figure 4 shows the predicted mixture relaxation modulus and shift factor mastercurves developed from the BBR test data. Figure 4a compares the mixtures’ relaxation modulus mastercurves of all four mixtures. The steeper slope over time indicates a higher relaxation rate of the sample and vice versa. Clearly, the 5% SBS modified mixture has the highest relaxation rate of the tested binders. The effective polymer network in the 5% SBS modified binder may contribute some additional relaxation ability to the binder. The 3% SBS modified mixture has slightly better relaxation rate than the base binder. On the other hand, the 1% SBS modified mixture has the lowest relaxation rate. In 1% SBS modified binder, an inadequate amount of polymer could not contribute to its elastic response improvement. Figure 4b shows the shift factors of different binders at different test temperatures.

Relaxation modulus and shift factor mastercurves of asphalt concrete mixtures: (a) relaxation modulus and (b) shift factor.
As mentioned above, this study performed IDT strength test at 0°C by applying a vertical constant crosshead displacement rate of 12.5 mm/min (0.208 mm/s). Equation 27 was used to determine the horizontal strain rate at the middle of the specimen. After obtaining the horizontal strain rate and the relaxation modulus, the IDT test data were used to determine the characteristic damage curve for each mixture. Figure 5 compares the IDT strength test results of all tested samples. Figure 5a compares the stress-strain relationships among the tested AC mixtures. It shows that the 5% SBS modified mixture has the highest IDT strength followed by the 3% SBS modified mixture. The increase in tensile strength of the 5% and 3% SBS modified mixture is attributed to the presence of effective polymer network in the binder. The base mixture has the lowest IDT strength. Though 1% SBS modified mixture has higher tensile strength, failure strain is the least among all mixtures. Figure 5b compares the characteristic damage curves for all mixtures. It reveals that damage resistance is higher for the 5% and 3% SBS modified mixtures. The 1% SBS modified mixture has the lowest resistance. After developing the damage curve of each AC mixture, they were used to calculate the tensile strength at thermal strain rate.

Indirect tension (IDT) test results at 0°C: (a) stress-strain plots and (b) characteristic damage curves.
Figure 6 compares the predicted thermal stress profile, failure stress, and critical cracking temperature with the measured data from TSRST. Figures 6, a–d, compare the predicted and measured thermal stresses for base, 1%, 3%, and 5% SBS modified mixtures, respectively. All these figures show that the predicted thermal stresses match well with the measured data. The R2 values for all cases are close to 1. Figure 6, e and f , compare the predicted critical cracking temperatures and failure strengths with the measured critical cracking temperatures and failure strengths, respectively. Figure 6e shows that the predicted critical cracking temperatures are close to the measured values for all cases (R2 = 0.85). Figure 6f shows that the two predicted strengths are close to the respective measured strengths and two predicted values are slightly higher than their counterparts. The R2 value for this case is 0.72, which also indicates a relatively good prediction. The percentage errors (absolute) between the measured and predicted values are shown in Table 3. The percentage errors are 3%–8% in failure temperature and 5%–14% in failure strength prediction. Given that all the percentage errors are below 15%, it can be said that a reasonable prediction can be made.
Percentage of Errors between Measured and Predicted Values

Comparison between measured and predicted data: (a) base mixture, (b) 1% SBS mixture, (c) 3% SBS mixture, (d) 5% SBS mixture, (e) measured versus predicted temperature, and (f) measured versus predicted strength.
Conclusion
This study presents a new mechanistic procedure to determine the critical cracking temperature of an AC mixture from BBR and IDT tests. The main reason for choosing these two tests is that they are widely used in the pavement industry for binder quality assurance and mix design purposes. This developed procedure uses TTSP, Prony series-based interconverted method, and Hirsch model to develop the mixture relaxation modulus mastercurve from the BBR data. The Hirsch model parameters need to be calibrated by comparing creep data from BBR and IDT creep tests performed at the same temperature. Next, it applies Boltzmann hereditary integral to calculate thermal stress over time caused by a constant rate of temperature drop. Before that, it applies the second-order heat equation to determine the actual temperature distribution in the sample. For the IDT strength data, the proposed procedure utilizes the VECD model to convert the IDT strength from test strain rate into thermal strain rate. Since it is not recommended to use any strain gauge for traditional IDT strength testing, this study derived an analytical equation based on the Hondros solution to compute the horizontal strain rate from the crosshead displacement rate. Finally, the critical cracking temperature is determined by coupling the thermal stress and IDT strength curves. Using the procedure presented in this paper, the critical cracking temperatures of four asphalt mixtures are predicted from BBR and IDT data. To validate the method, their actual critical cracking temperatures were measured using TSRST performed in the laboratory. The predicted critical cracking temperatures are found to be very close to the laboratory measured values. The developed procedure will have substantial practical and technical importance in predicting the critical cracking temperature of an AC mixture with reasonable accuracy.
In the future, a full-scale parametric study will be pursued to evaluate the sensitivity of each model parameter for predicting the critical cracking temperature. In addition, the results from the proposed method will be compared with the results from other methods, like the dissipated pseudo-strain energy-based failure criterion developed by Keshavarzi and Kim ( 35 ).
Footnotes
Acknowledgements
The authors would like to express their sincere gratitude and appreciations to the Project Technical Panel Members, Project Advocate, Project Sponsor, and Project Manager at the NMDOT.
Author Contributions
The authors confirm the contribution to the paper as follows: study conception and design: M. A. Hasan, R. A. Tarefder; data collection: M. A. Hasan; analysis and interpretation of results: M. A. Hasan, R. A. Tarefder; draft manuscript preparation: M. A. Hasan, R. A. Tarefder. All authors reviewed the results and approved the final version of the manuscript.
Declaration of Conflicting Interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This study was funded by the New Mexico Department of Transportation (NMDOT).
