Advances in Permanent Deformation Modeling of Asphalt Concrete—A Review

Permanent deformation is one of the dominant asphalt concrete damages. Significant progress has been made to realistically predict the damage. In the last decade, the mechanistic approach has been the focus of research, and the fundamental theories of viscoelasticity, viscoplasticity, continuum mechanics, and micromechanics are applied to develop the material laws (constitutive equations). This paper reviews the advancement of permanent deformation models including analogical, microstructural, and continuum-based methods. Pavement analysis using the nonlinear damage approach (PANDA) is the most comprehensive and theoretically sound approach that is available in the literature. The model coupled different damages and other phenomena (such as cracking, moisture, and phenomena such as healing, aging, etc.). The anisotropic microstructure approach can be incorporated into the PANDA approach for a more realistic prediction. Moreover, the interaction of fatigue and permanent deformation is the gap that is lacking in the literature. The mechanistic approaches have the capacity to couple these damages for unified asphalt concrete damage prediction.


Introduction
A flexible pavement, the longest continuous structure, comprises an asphalt concrete layer supported by unbound compacted layers (aggregate bases and subgrade). Asphaltic materials have been used for roadway construction since the end of the nineteenth century [1]. Asphalt concrete (also called bituminous mixture or hot mix asphalt) is a complex heterogeneous and three-phase material (aggregate matrix, mastic, and the air void). Such a material's performance depends on the mixture composition, proportion, mechanical properties, and environmental conditions. It is characterized as a viscoelastic, viscoplastic, and time-and temperature-dependent material. Due to external loads and environmental factors, different distresses or damages occur in the asphalt layer. Permanent deformation or rutting is one of the primary distresses, making pavements rough and unsafe for driving, causing hydroplaning, etc. The rutting distress was noticed as a primary asphalt performance criterion [2]. The asphalt concrete's susceptibility to permanent deformation is linked to material attributes and climatic and loading factors [3][4][5]. Material-related factors include excessive asphalt content, fine aggregate, high natural sand percentage, rounded aggregate particles, the moisture content in the mix, or granular materials and soils. From the asphalt concrete constituent properties, the chemistry of asphalt binder is the component (if not only) that makes bituminous mixtures a complex rate-dependent, nonlinear material. This nonlinear (viscous) behavior of the binder makes the permanent deformation evolution of asphalt mixtures a nonlinear mechanism [6]. For this reason, the rheological characteristics of binder (shear modulus, G*, and viscosity, η) are used to classify deformability properties of mixtures. The Strategic Highway Research Program (SHRP) rheological parameter, G*/sin δ (δ-the phase angle), is the widely used changeably. The term "Rutting" is used to describe the pavement surface r to the vertical depression along the wheel path caused by the permanent d wear in the asphalt layer. Permanent deformation is the accumulation of strain due to densification, shear deformation, and crack growth in aspha this paper, the term permanent deformation is used to refer to the plastic an strains (irrecoverable deformation). Therefore, rutting (RD, mm) is express where h is the ith layer's thickness, and ε , is the viscoplastic/permanen ith layer.

Objective and Scope
The aim of this review is to provide state-of-the art information on the of permanent deformation modeling for asphalt concrete. In the review, th nent deformation modeling theories, methods, models, and calibration tests in detail. The focus of the paper is studying the advancement of constitu approaches (analogical, microstructural, and continuum-based) and synthe ity of the approaches, merits, and limitations. The robustness of permanen models to account for/couple simultaneous damages such as moisture, fat etc., were also explored. The literature studied in this paper was collected u (strings) such as viscoelastic, viscoplastic, viscodamage, permanent deform continuum damage, microstructure or micromechanics, mechanistic meth covery, etc. The organization of the paper is depicted in the flow chart in Fi

Analytical Models
The evolution of the irrecoverable deformation of asphaltic materia loading is described by three distinct stages, as shown in Figure 2a

Analytical Models
The evolution of the irrecoverable deformation of asphaltic material due to cyclic loading is described by three distinct stages, as shown in Figure 2a. The primary zone is described by the rapid accumulation of permanent deformation at a decreasing strain rate. In the secondary zone is the constant rate of permanent deformation with strain hardening, and the tertiary stage is characterized by the increasing rate of deformation accumulate and crack formation. The flow time (FT) or the flow number (FN) is defined as the time or number of loading cycles when shear deformation under constant volume commences. Several researchers have verified that the asphalt concrete's deformation evolution showed all the three phases [28][29][30][31][32]. This three-stage deformation property is also regarded as asphalt concrete material property. The most common and simulative test used to characterize permanent deformation is the triaxial creep-recovery test. An example in Figure 2b shows that the creep-recovery deformation is dependent on confining stress. Confinement increases the friction between aggregates and increases the resistance to deformation. Although there is no conclusive research, the in situ confining pressure of asphalt concrete is approximated to be between 100 kPa and 225 kPa.
OR PEER REVIEW 4 of 29 hardening, and the tertiary stage is characterized by the increasing rate of deformation accumulate and crack formation. The flow time (FT) or the flow number (FN) is defined as the time or number of loading cycles when shear deformation under constant volume commences. Several researchers have verified that the asphalt concrete's deformation evolution showed all the three phases [28][29][30][31][32]. This three-stage deformation property is also regarded as asphalt concrete material property. The most common and simulative test used to characterize permanent deformation is the triaxial creep-recovery test. An example in Figure 2b shows that the creep-recovery deformation is dependent on confining stress. Confinement increases the friction between aggregates and increases the resistance to deformation. Although there is no conclusive research, the in situ confining pressure of asphalt concrete is approximated to be between 100 kPa and 225 kPa. Several analytical models have been proposed over several decades to predict the permanent deformation of asphalt concrete. Sousa et al. [33] gave a summary of selected analytical models. Some of the commonly used permanent deformation models and the corresponding calibration tests are summarized in Table 1. These analytical models are used to predict rutting from typical laboratory experimental data (axial stress-strain test data). The models in Table 1 can be classified as empirical and mechanistic-empirical, which are calibrated using simulative laboratory or field data. The models give the macroscopic responses of the measured data and hardly relate the fundamental material properties. Other models are regression equations to fit rutting data and lack the explicit physical or material property for the modeling parameters.

Calibration Tests
The most common laboratory test protocols used cylindrical specimens of dimensions (diameter by height) 100 mm by 150 mm for creep-recovery, or 150 mm by 50-70 mm for creep. The shear strain from simple shear tests is also used to model permanent deformation [16,34]. The coupling of shear and axial strain components in permanent deformation modeling has not been performed yet in previous studies [35]. The indirect tensile test is also used for permanent deformation with different specimen dimensions of 150 or 100 mm diameter by 50 to 70 mm thickness [36]. Moreover, the wheel tracking test is used to simulate permanent deformation [37][38][39]. Several analytical models have been proposed over several decades to predict the permanent deformation of asphalt concrete. Sousa et al. [33] gave a summary of selected analytical models. Some of the commonly used permanent deformation models and the corresponding calibration tests are summarized in Table 1. These analytical models are used to predict rutting from typical laboratory experimental data (axial stress-strain test data). The models in Table 1 can be classified as empirical and mechanistic-empirical, which are calibrated using simulative laboratory or field data. The models give the macroscopic responses of the measured data and hardly relate the fundamental material properties. Other models are regression equations to fit rutting data and lack the explicit physical or material property for the modeling parameters.

Calibration Tests
The most common laboratory test protocols used cylindrical specimens of dimensions (diameter by height) 100 mm by 150 mm for creep-recovery, or 150 mm by 50-70 mm for creep. The shear strain from simple shear tests is also used to model permanent deformation [16,34]. The coupling of shear and axial strain components in permanent deformation modeling has not been performed yet in previous studies [35]. The indirect tensile test is also used for permanent deformation with different specimen dimensions of 150 or 100 mm diameter by 50 to 70 mm thickness [36]. Moreover, the wheel tracking test is used to simulate permanent deformation [37][38][39]. Table 1. Some selected permanent deformation analytical models.

Model (Equation)
Variables Description Reference Calibration Test Accurate for secondary stages, small stress/strain deformation Creep, creep-recovery Most widely used analytical model for permanent deformation in all three creep stages; the first part is a power function (for low stresses). and the second part is for high stresses (tertiary stage). (σ VM -maximum stress, σ VL -plastic failure threshold; σ 1 , σ 3 -axial and lateral stress, E * -stiffness) Francken [40] Repeated triaxial compression Analytical Power Model based on dissipated energy rate; A is a function of resilient modulus and applied stress Khedr Safwan [41] Multiple step dynamic test Analytical model for three creep stages, mainly developed for unbound materials (δ 1 , δ 3 -scale primary and tertiary strain; δ 2 , δ 4 -rate parameters) Wilshire and Evans [43] Creep tests  Wheel tracking, Uniaxial cyclic compression As shown in Figure 3, different models have different accuracy for the same permanent deformation data. It is evident that the fitting accuracy is variable especially in the primary and tertiary stages of deformation. The incremental or Choi, Tseng-Lytton, and Francken models showed close predictions of the measured data. The Francken model is the most widely used for permanent deformation modeling [30].
Each model has a different accuracy of fitting all the creep stages of permanent deformation. Some of the models presented above have clear limitations such as the implicit empiricism, being unable to model load history and hardening-relaxation behavior, a lack of capacity to couple other simultaneous damages, etc. As discussed in the next sections, mechanistic permanent deformation prediction methods are aimed to resolve the limitations of empirical/mechanistic-empirical models. The latest mechanistic models have relied on rigorous material models, which are based on fundamental theories of mechanics, stress-strain relationships, and environmental factors. Francken models showed close predictions of the measured data. The Francken model is the most widely used for permanent deformation modeling [30].
(a) (b) Each model has a different accuracy of fitting all the creep stages of permanent deformation. Some of the models presented above have clear limitations such as the implicit empiricism, being unable to model load history and hardening-relaxation behavior, a lack of capacity to couple other simultaneous damages, etc. As discussed in the next sections, mechanistic permanent deformation prediction methods are aimed to resolve the limitations of empirical/mechanistic-empirical models. The latest mechanistic models have relied on rigorous material models, which are based on fundamental theories of mechanics, stress-strain relationships, and environmental factors.

Stress-Strain Response
The stress-, time-, and temperature-dependent viscoelastic, viscoplastic, and viscodamage properties of asphalt concrete material offered considerable challenges to accurately model the response under variable loading conditions. The response of asphalt concrete is stress path-dependent [26]. It is also a rate-and history-dependent material [47,48]. The linearity limits of asphalt concrete are 150 and 100 micro-strain in compression and tension, respectively [49], but others suggest 122 micro-strain as a limit for linear response [50]. Permanent deformation damage is induced on the asphalt beyond this strain limit at high temperatures. Moreover, the stress-strain evolution is highly dependent on the stress and strain levels, number of loading cycles (loading time), and temperature range. In a typical creep-recovery test, the stress-strain hysteresis loops evolve nonlinearly as shown in Figure 4a. The loop has the recovery and non-recovery (permanent deformation) parts. In each creep-recovery cycle, the deviatoric stress causes a non-recoverable strain and creates an open stress-strain hysteresis loop. Figure 4b shows the stress-strain responses of a constant rate compressive strength test (without recovery time). Asphalt concrete undergoes creep deformation during the load phase and a delayed recovery upon the load removal during the rest period. Traditionally, the additive decomposition of strain is applied to separate the permanent strain and recoverable strain from the total strain.

Stress-Strain Response
The stress-, time-, and temperature-dependent viscoelastic, viscoplastic, and viscodamage properties of asphalt concrete material offered considerable challenges to accurately model the response under variable loading conditions. The response of asphalt concrete is stress path-dependent [26]. It is also a rate-and history-dependent material [47,48]. The linearity limits of asphalt concrete are 150 and 100 micro-strain in compression and tension, respectively [49], but others suggest 122 micro-strain as a limit for linear response [50]. Permanent deformation damage is induced on the asphalt beyond this strain limit at high temperatures. Moreover, the stress-strain evolution is highly dependent on the stress and strain levels, number of loading cycles (loading time), and temperature range. In a typical creep-recovery test, the stress-strain hysteresis loops evolve nonlinearly as shown in Figure 4a. The loop has the recovery and non-recovery (permanent deformation) parts. In each creep-recovery cycle, the deviatoric stress causes a non-recoverable strain and creates an open stress-strain hysteresis loop. Figure 4b shows the stress-strain responses of a constant rate compressive strength test (without recovery time). Asphalt concrete undergoes creep deformation during the load phase and a delayed recovery upon the load removal during the rest period. Traditionally, the additive decomposition of strain is applied to separate the permanent strain and recoverable strain from the total strain.  The schematic Figure 5 shows the strain components of a single pulse creep-recovery loading. The strain components can be separated into four strain components [51,52] (elastic ε , viscoelastic ε , plastic ε , and viscoplastic ε ). The elastic (time-independent) and viscoelastic (time-dependent) are recoverable, while plastic (time-independent) and viscoplastic (time-dependent) strains are non-recoverable parts of total strain. The total strain is expressed as follows. The schematic Figure 5 shows the strain components of a single pulse creep-recovery loading. The strain components can be separated into four strain components [51,52] (elastic ε e , viscoelastic ε ve , plastic ε p , and viscoplastic ε vp ). The elastic (time-independent) and viscoelastic (time-dependent) are recoverable, while plastic (time-independent) and viscoplastic (time-dependent) strains are non-recoverable parts of total strain. The total strain is expressed as follows. The schematic Figure 5 shows the strain components of a single pulse creep-recovery loading. The strain components can be separated into four strain components [51,52] (elastic ε , viscoelastic ε , plastic ε , and viscoplastic ε ). The elastic (time-independent) and viscoelastic (time-dependent) are recoverable, while plastic (time-independent) and viscoplastic (time-dependent) strains are non-recoverable parts of total strain. The total strain is expressed as follows.
The permanent strain (ε = ε + ε ) is obtained by subtracting the viscoelastic strain (ε = ε + ε ) components from the total strain. However, the strain decomposition approach is questioned when the hardening-relaxation mechanism of asphalt concrete is considered. In a loading-unloading cycle, the recoverable strain is dependent on the rest period. For tests with short rest periods, the computed viscoplastic strain can be overestimated [53]. The limitation of strain decomposition is from the inherent interaction between the viscoelastic and viscoplastic strain (the viscoelastic response is also a function of viscoplastic deformation history).

Constitutive Models
The computational modeling of asphalt concrete poses difficulties mainly due to the material nonlinearity, complexity to characterize under repeated and moving loads, and variable environmental conditions (temperature, moisture, etc.) [35,54]. The constitutive equations for the linear viscoelastic strain (ε ) in undamaged conditions is defined by the Boltzmann superposition principle. The permanent strain (ε vp = ε p + ε vp ) is obtained by subtracting the viscoelastic strain (ε ve = ε e + ε ve ) components from the total strain. However, the strain decomposition approach is questioned when the hardeningrelaxation mechanism of asphalt concrete is considered. In a loading-unloading cycle, the recoverable strain is dependent on the rest period. For tests with short rest periods, the computed viscoplastic strain can be overestimated [53]. The limitation of strain decomposition is from the inherent interaction between the viscoelastic and viscoplastic strain (the viscoelastic response is also a function of viscoplastic deformation history).

Constitutive Models
The computational modeling of asphalt concrete poses difficulties mainly due to the material nonlinearity, complexity to characterize under repeated and moving loads, and variable environmental conditions (temperature, moisture, etc.) [35,54]. The constitutive equations for the linear viscoelastic strain (ε ve ) in undamaged conditions is defined by the Boltzmann superposition principle.
where t is time; E(t) and D(t) are relaxation and creep compliance moduli, respectively; and τ is the integration variable. Prony series forms of the creep compliance and relaxation moduli are expressed as follows.
where N and M are the total numbers of Prony terms; D o , D i , and τ i are creep compliance model coefficients; and E ∞ , E i , and ρ i are relaxation model coefficients. Often, D(t) is obtained from E(t) data via the interconversion technique [55,56] . Once the creep compliance function is defined, the viscoelastic strain can be determined, and the viscoplastic strain is calculated from the total strain by the additive decomposition technique. A sufficient rest period is necessary to completely remove the delayed recovery strain from viscoplastic strain evolution. Cao and Kim [53] showed that 99% of viscoelastic strain is recovered for the 0.4 s pulse period and 100 s rest duration in the first cycle and 98% in cycle number 10 within 60 s. Therefore, to obtain a true viscoplastic strain, about 100 s rest period is required [45]. Experimental observations showed that the total deformation in each cycle decreases due to hardening as the number of load cycles increases. Since the viscoelastic strain is obtained from a separate dynamic modulus test (i.e., a constant viscoelastic strain), subtracting a constant cyclic viscoelastic strain from a decreasing total strain can result in a negative viscoplastic strain. That means a decreasing viscoelastic deformation model should be proposed. It is the microstructural change due to viscoplastic deformation that causes a change in viscoelastic deformation. This interaction between viscoelastic and viscoplastic deformation is referred to as viscoelastic-viscoplastic coupling, according to [53]. There is no available literature that couples the two deformations.

Analogical Models
The family of different analogical models has been used to model the viscoelasticviscoplastic response of time-dependent materials [57][58][59][60]. The common classic mechanical models are spring, dashpot, slip device, pot, parabolic elements, etc., and the combination of these elements to analogs. The mechanical elements are advantageous to visualizing the stress and strain responses using the analogs. The Maxwell model (for the viscoelastic model), Kelvin model (for creep response), and the Burgers model are used to model the viscoelastic and viscoplastic strain (Figure 6a,b). The governing differential equations (viscoelastic constitutive equations) were developed from a number of springs and dashpots arranged in series and parallel. In the generalized Maxwell model, the same strain is shared across all elements, and the stress is additive, while in the generalized Burgers model the strains are additive, and the stress is the same for each element. It can be noted here that the generalized Burgers model shares the same framework as classical viscoplasticity models and allows nonlinearities based on stress to be accommodated more easily [57]. The viscoelastic and viscoplastic components can be calculated using the hereditary integral formulation as follows.
D and D are viscoelastic and viscoplastic creep compliance. The hereditary integrals in Equations (7) and (8) are different from the one in Equation (3). A formulation D ve and D vp are viscoelastic and viscoplastic creep compliance. The hereditary integrals in Equations (7) and (8) are different from the one in Equation (3). A formulation based on stress and the rate of compliance rather than a formulation based on the rate of stress and compliance is advantageous to avoid problems due to the sudden application of a stress in which the rate of stress can be extremely high (e.g., in a creep test). The first derivatives of the viscoelastic and viscoplastic creep compliance for a generalized Burg- where λ i is the viscosity of the ith Voigt element; E i is the modulus of elasticity; and λ ∞ is the viscosity of viscoplastic element. A power function for creep compliance is also used for small stress cases. Moreover, the model in Figure 7c is a modified Burgers model with a plasticity element for asphalt mixture [60,61]. The additional elastoplastic network composed of the spring and slider in parallel is used. The limit stress in the slider modeling plasticity is denoted by σ o . The authors also extended the fractional rheological model for nonlinear elastic, nonlinear viscous, and plastic properties and formulated a differential equation to characterize the viscoelastic-plastic response of asphalt concrete. Other similar analogical models such as the 2S2P1D (two springs, two parabolic. and one dashpot elements) in Figure 7a and DBN (Di Benedetto-Neifar) in Figure 7b are also frequently used to predict linear viscoelastic and creep responses (for a small number of cycles) for binders and bituminous mixtures [58,59,62]. The DBN model is a special case of the Kelvin-Voigt model where the DBN model has an elastoplastic (EP) element instead of an elastic element only. From these models, the governing differential equation for creep compliance or stress and strain functions is derived to predict the permanent deformation response of asphalt concrete.
One can note that the parabolic elements (k, h) in 2S2P1D model, the slider element in (Figure 6c), and the spring-pot element in the fractional model ( Figure 6d) have similar arrangements. The parabolic creep elements in 2S2P1D (Figure 7a), the elastic-plastic (EP) elements in DNB (Figure 7b) and the spring-pot elements in fractional model (Figure 6d) have also similar functions to model the elastoplastic response of asphalt concrete. The slip device shown in Figure 7c is placed in parallel with linear viscoelastic (LVE) element which has similar property as the fractional model in Figure 6c. The slip device functions as an irreversible deformation, which means that whenever the LVE device relaxes during unloading, the slip device locks, thereby disallowing strain recovery. Based on this phenomenology, the viscoelastic integrals are proposed for hardening and permanent deformation, considering that viscoelastic deformation, viscoplastic deformation, and hardening function are history dependent. Subramanian et al. [48] proposed a viscoelastic-like viscoplastic constitutive model for the permanent deformation of asphalt concrete. The proposed model takes the following form in Macaulay brackets. From these models, the governing differential equation for creep compliance or stress and strain functions is derived to predict the permanent deformation response of asphalt concrete.

dσ (τ)
One can note that the parabolic elements (k, h) in 2S2P1D model, the slider element in (Figure 6c), and the spring-pot element in the fractional model ( Figure 6d) have similar arrangements. The parabolic creep elements in 2S2P1D (Figure 7a), the elastic-plastic (EP) elements in DNB (Figure 7b) and the spring-pot elements in fractional model (Figure 6d) have also similar functions to model the elastoplastic response of asphalt concrete. The slip device shown in Figure 7c is placed in parallel with linear viscoelastic (LVE) element which has similar property as the fractional model in Figure 6c. The slip device functions as an irreversible deformation, which means that whenever the LVE device relaxes during unloading, the slip device locks, thereby disallowing strain recovery. Based on this phenomenology, the viscoelastic integrals are proposed for hardening and permanent deformation, considering that viscoelastic deformation, viscoplastic deformation, and hardening function are history dependent. Subramanian et al. [48] proposed a viscoelastic-like viscoplastic constitutive model for the permanent deformation of asphalt concrete. The proposed model takes the following form in Macaulay brackets.
where H(t) is the material hardening variable, H o denotes the initial hardening state, σ H (t) and σ vp (t) are functions of deviatoric stress for hardening and viscoplastic deformation calculation, and D 1 (t) and D 2 (t) are compliance functions. The two stress terms are approximated by power functions, as follows.
where σ d (t) is the deviatoric stress history; and H 1 , q 1 , G 1 , G 2 , p 1 , p 2 , and α are parameters. This model considers only hardening during loading pulses and ignored the softening mechanism. Based on the model in [48], Cao and Kim [53] proposed a viscoplastic model using "internal stress" as the hardening variable. They hypothesized that as soon as the deformation of the LVE device is constrained (in Figure 7c), the internal stress inside the device starts to develop due to stress relaxation. As illustrated in Figure 8, the internal stress decreases with time due to LVE device relaxation before the next load cycle but has a direction opposite to the external load. Once the applied stress rises above the level of the concurrent resisting internal stress in the LVE device, the slip device is unlocked and becomes frictionless, allowing the overall deformation of the mechanical analog to increase in a viscoelastic fashion. The proposed viscoelastic-type viscoplastic model takes the following form [53].
level of the concurrent resisting internal stress in the LVE device, the slip device is unlocked and becomes frictionless, allowing the overall deformation of the mechanical analog to increase in a viscoelastic fashion. The proposed viscoelastic-type viscoplastic model takes the following form [53]. Adapted with permission from [53]. 2016, Mechanics of Materials, Elsevier.
This model introduces the coupling of viscoelastic and viscoplastic responses using the internal stress as a hardening variable and viscoelastic-like hardening-relaxation spectrum.

Continuum Based Models
Continuum mechanics is a standalone and widely applied theory for damage formulation. The concept of continuum damage mechanics (CDM) was pioneered by Kachanov This model introduces the coupling of viscoelastic and viscoplastic responses using the internal stress as a hardening variable and viscoelastic-like hardening-relaxation spectrum.

Damage Density
Continuum mechanics is a standalone and widely applied theory for damage formulation. The concept of continuum damage mechanics (CDM) was pioneered by Kachanov [63], who introduced a scalar measure called damage variable or damage density φ, which is defined as follows.
where A-real (intact) area, A d -damaged area, A d − A is the area of micro-damage, φ = 0 means the initial state, and φ = 1 mean complete rupture. Based on the damage density function and effective area, the effective stress concept in CDM is defined as follows.
where σ ij is the effective stress tensor in an undamaged configuration; and σ ij is the nominal Cauchy tensor in damage configuration. For more accurate modeling, the damage evolution is modified as follows [64].
The classic Kachanov-Robotnov damage models [63,65] were extensively used for creep damage (φ c ) modeling for different materials.
where Γ ϕ o is reference damage viscosity, Y 0 is reference damage force, Y is damage driving force in effective configuration, G(T) is temperature coupling term, and k and q are constants.

Viscoplasticity
The classic Perzyna viscoplastic hardening rule [66] assumes a constant hardening variable for a cyclic creep-recovery load, defined as where Γ vp is the viscoplastic fluidity parameter such that 1/Γ vp is a measure of viscoplastic viscosity; and ∂σ is a measure of the direction of viscoplastic strain. The classical hardening assumes that the viscoplastic strain rate decreases progressively with an increase in loading time. However, the hardening function is not constant due to the hardeningrelaxation behavior [23]. Different researchers pointed out that the Perzyna-type rate models have limitations such as that the model cannot capture the load history effect, the relaxation or softening behavior during the rest period is ignored, and it assumes a constant hardening parameter [23,47,48]. As illustrated in Figure 9, the viscoplastic strain rate is no longer a decreasing function. The quantity q vp is the hardening-relaxation internal state variable that memorizes the maximum experienced viscoplastic strain for which the hardening recovery has occurred. ening-relaxation behavior [23]. Different researchers pointed out that the Perzyna-type rate models have limitations such as that the model cannot capture the load history effect, the relaxation or softening behavior during the rest period is ignored, and it assumes a constant hardening parameter [23,47,48]. As illustrated in Figure 9, the viscoplastic strain rate is no longer a decreasing function. The quantity q vp is the hardening-relaxation internal state variable that memorizes the maximum experienced viscoplastic strain for which the hardening recovery has occurred.  Pavement analysis using the nonlinear damage approach (PANDA) is the latest generation of mechanistic pavement design approach [67]. The PANDA is a mechanistic-based pavement analysis method that is founded on three classic theories: (1) Schapery's [68] nonlinear viscoelasticity, (2) Perzyna's [66] viscoplasticity, and (3) Darabi's [64] viscodamage constitutive relationship. Based on the three constitutive equations, the PANDA approach has an unlimited capacity to couple different damage mechanisms of pavement structures. The approach coupled the temperature, rate, and time-dependent viscoelastic and viscoplastic models to predict the permanent deformation of asphalt concrete. For example, healing, aging, hardening-relaxation, moisture-induced damage, and other behaviors are conveniently incorporated into the PANDA approach [69][70][71][72][73]. In Figure 10, the development of the PANDA model and the constitutive equations are summarized.
(1) First, the linear and nonlinear viscoelastic variables are obtained from dynamic modulus (for linear viscoelastic) and creep-recovery tests (nonlinear viscoelastic). The nonlinear viscoelastic strain is formulated using the well-known Schapery's viscoelastic constitutive equation [68]. (2) Secondly, the viscoelastic strain is deducted from the total strain to extract the viscoplastic strain from the same creep-recovery test data using the strain decomposition principle. Then, the classic Perzyna's viscoplasticity [66] theory is adopted to predict the viscoplastic strain evolution. The Drucker-Prager yield surface function is often used [23,74].
(3) The third foundation of the PANDA constitutive model is the viscodamage model [64] using the continuum damage mechanics (CDM) theory. The effective strain is used in effective configuration. The viscodamage model mainly predicts the permanent deformation in the tertiary creep.
Therefore, the PANDA model encompasses nonlinear viscoelastic, viscoplastic (hardening), and viscodamage responses. The thermo-piezo-rheological viscoelastic properties [75] coupled with the viscoplastic yield criteria (the Drucker-Prager yield surface) is an integral part of the PANDA method [76]. Several material parameters need to be optimized to calibrate the viscoelastic-viscoplastic-viscodamage, the hardening-relaxation, moisture damage, and healing responses of asphalt concrete. The parameters and their physical meaning are presented in Table 2 Despite the unlimited capacity of the PANDA approach, calibrating the large number of model parameters is a laborious task. Hence, a systematic procedure is followed to obtain material parameters with a smaller number of tests. Once robust mechanistic constitutive models are developed and calibrated, the numerical implementation (finite element modeling) is performed using the user-defined material (UMAT) tool to define material laws in commercial software such as ABAQUS. Finite element modeling (FEM) is conducted at the desired modeling space (2D or 3D), and realistic tire-pavement contact [77], traffic, and full pavement structure (asphalt, base, subbase, and subgrade), etc., can be constructed [67,76,[78][79][80]. Although the PANDA approach is still at a research stage, it is evident that advantages as well as limitations can be listed. In Table 3, some of the merits and limitations of the PANDA are described. The calibration of PANDA models used uniaxial creep, uniaxial constant stress creep-recovery, the crosshead strain rate test, and multiple stress creep tests, etc. The influence of confining pressure on the linear viscoelastic as well as nonlinear viscoelastic responses is significant [75,81]. Most studies that tried to calibrate the PANDA models used uniaxial test data, or some used a single confining pressure. Therefore, the PANDA model encompasses nonlinear viscoelastic, viscoplastic (hardening), and viscodamage responses. The thermo-piezo-rheological viscoelastic properties [75] coupled with the viscoplastic yield criteria (the Drucker-Prager yield surface) is an integral part of the PANDA method [76]. Several material parameters need to be optimized to calibrate the viscoelastic-viscoplastic-viscodamage, the hardening-relaxation, moisture damage, and healing responses of asphalt concrete. The parameters and their physical meaning are presented in Table 2 Despite the unlimited capacity of the PANDA approach, calibrating the large number of model parameters is a laborious task.

Microstructural Based Models
The micromechanics approach is probably the best way to account the effects of individual mixture constituents and their interactions and the anisotropy of heterogenous asphalt mixture. The microstructure change in asphalt concrete is mainly due to the friction between the aggregates and interlocking bond breakage. This mechanism is responsible for the accumulation of permanent deformation rather than the recoverable part of the deformation. The continuous increase in the resistance of the material due to the permanent microstructure rearrangement is physically related to the strain hardening. The hardening parameter reflects the combined effect of the cohesion of asphalt binders, the adhesion properties between binder and aggregate, and the frictional properties of the aggregate structure. The fabric of granular media refers to the size, shape, and arrangement of the solid particles and the associated voids. The scalar quantity, like void ratio, is not capable of characterizing the directional nature of fabric and describing the state of packing of granular materials [85]. The microstructure approach is necessary to consider the directional nature of granular fabric. The approach is capable of modeling nonlinearities such as heterogeneity, aggregate distribution, anisotropy, crack, and air void. The micromechanics coupled with the continuum damage approach is a powerful way to model the permanent deformation of granular asphalt concrete with an appropriate representative volume element (RVE) [19][20][21][22]. The fundamental element of the granular microstructure is a directional unit vector m and the vector magnitude ∆ (Figure 11).
where M is the total number of objects analyzed in an image; and θ k is the orientation of the unit vector n and ranges between −90 and +90 • . Theoretically, the value of ∆ is between zero and unity; ∆ = 0 indicates that objects are completely randomly distributed, which is analogous to isotropic materials; and ∆ = 1 indicates that objects are entirely oriented in one direction, which is analogous to perfectly transverse anisotropic materials. Oda and Nakayama [86] introduced a symmetrical microstructure tensor F ij that gives a measure of the two-dimensional anisotropy produced by the preferred orientation of non-spherical particles. 23) state of packing of granular materials [85]. The microstructure approach is necessary to consider the directional nature of granular fabric. The approach is capable of modeling nonlinearities such as heterogeneity, aggregate distribution, anisotropy, crack, and air void. The micromechanics coupled with the continuum damage approach is a powerful way to model the permanent deformation of granular asphalt concrete with an appropriate representative volume element (RVE) [19][20][21][22]. The fundamental element of the granular microstructure is a directional unit vector m and the vector magnitude ∆ (Figure 11). Figure 11. Two-dimensional particles orientation for micromechanics modeling (orientation, and vector, m).
where M is the total number of objects analyzed in an image; and is the orientation of the unit vector n and ranges between −90 and +90°. Theoretically, the value of ∆ is between zero and unity; ∆ = 0 indicates that objects are completely randomly distributed, which is analogous to isotropic materials; and ∆ = 1 indicates that objects are entirely oriented in one direction, which is analogous to perfectly transverse anisotropic materials. Oda and Nakayama [86] introduced a symmetrical microstructure tensor that gives a measure of the two-dimensional anisotropy produced by the preferred orientation of nonspherical particles. The microstructure tensor is incorporated in the Drucker-Prager yield function by modifying the first and second invariants. The classical Perzyna's viscoplastic and continuum damage models are then used to formulate the microstructure-based viscoplastic strain model to characterize permanent deformation of asphalt concrete. The detail derivation of the model can be found in the study [19,20,86]. The final form of the viscoplastic model in a triaxial compression test (axial strain) based on micro-structural anisotropy takes the following form. For simplicity the anisotropy value can be taken as constant (the initial value). For isotropic elements, F 11 = (−4∆)/3(3 + ∆) and The microstructure tensor is incorporated in the Drucker-Prager yield function by modifying the first and second invariants. The classical Perzyna's viscoplastic and continuum damage models are then used to formulate the microstructure-based viscoplastic strain model to characterize permanent deformation of asphalt concrete. The detail derivation of the model can be found in the study [19,20,86]. The final form of the viscoplastic model in a triaxial compression test (axial strain) based on micro-structural anisotropy takes the following form.
ε vp 11 (24) where , and λ and µ are anisotropy coefficients that reflect the effect of the aggregate anisotropic distribution on the confining and shear stresses, respectively. The viscoplastic model considers phenomena including the effect of the binder fluidity (Γ); confinement and aggregate friction (α); aggregate interlocking and dilation (β); binder cohesion and its adhesion to the aggregates (κ); anisotropy of aggregate distribution (∆); and microstructure damage (ξ). As shown in Figure 12, the anisotropy parameter significantly contributes to compressive viscoplastic behavior and has little effect or no effect on the tensile and linear viscoelastic property [22]. In a conventional continuum method (isotropic mode), phenomenologically motivated, microstructuralbased viscoplastic models have been developed, for example in [23,25] and others. effect of the binder fluidity (Γ); confinement and aggregate friction ( ); aggregate interlocking and dilation ( ); binder cohesion and its adhesion to the aggregates ( ); anisotropy of aggregate distribution (∆); and microstructure damage ( ). As shown in Figure 12, the anisotropy parameter significantly contributes to compressive viscoplastic behavior and has little effect or no effect on the tensile and linear viscoelastic property [22]. In a conventional continuum method (isotropic mode), phenomenologically motivated, microstructural-based viscoplastic models have been developed, for example in [23,25] and others. Microstructure modeling is directly integrated with the utilization of digital technologies to capture the particles arrangement, deformation, etc., in the granular materials packing. Digital image correlations, X-ray chromatography, and other tools were used for asphalt concrete. Coleri et al. [88] used the X-ray computed tomography (CT) method to model asphalt concrete rutting in a full-scale heavy vehicle simulation (HVS) test site. The object segmentation and processing procedure using a digital camera is shown in Figure  13, and an example of the discretization phases is shown in Figure 14. Using the digital technology, parts of a heterogenous mixture can be easily dissected, modeled, and analyzed into the FEM and other post processes. Microstructure modeling is directly integrated with the utilization of digital technologies to capture the particles arrangement, deformation, etc., in the granular materials packing. Digital image correlations, X-ray chromatography, and other tools were used for asphalt concrete. Coleri et al. [88] used the X-ray computed tomography (CT) method to model asphalt concrete rutting in a full-scale heavy vehicle simulation (HVS) test site.
The object segmentation and processing procedure using a digital camera is shown in Figure 13, and an example of the discretization phases is shown in Figure 14. Using the digital technology, parts of a heterogenous mixture can be easily dissected, modeled, and analyzed into the FEM and other post processes. Microstructure modeling is directly integrated with the utilization of digital technologies to capture the particles arrangement, deformation, etc., in the granular materials packing. Digital image correlations, X-ray chromatography, and other tools were used for asphalt concrete. Coleri et al. [88] used the X-ray computed tomography (CT) method to model asphalt concrete rutting in a full-scale heavy vehicle simulation (HVS) test site. The object segmentation and processing procedure using a digital camera is shown in Figure  13, and an example of the discretization phases is shown in Figure 14. Using the digital technology, parts of a heterogenous mixture can be easily dissected, modeled, and analyzed into the FEM and other post processes.  In general, the microstructural based micromechanics modeling is ideal to consider the effect of individual constituent particles in viscoplastic evolution [89]. However, this approach needs image analyzing tools, discretization, and finite element modeling. Several material behaviors such as hardening-relaxation, healing, stress history, viscodamage, etc., need to be integrated for accurate formulation. Limited literature can be found in this area. With the advancement of morphological image analyzing tools, the In general, the microstructural based micromechanics modeling is ideal to consider the effect of individual constituent particles in viscoplastic evolution [89]. However, this approach needs image analyzing tools, discretization, and finite element modeling. Several material behaviors such as hardening-relaxation, healing, stress history, viscodamage, etc., need to be integrated for accurate formulation. Limited literature can be found in this area. With the advancement of morphological image analyzing tools, the microstructural approach will be the next active research area to characterize permanent deformation.

Finite Element Simulation
Even if pavement is considered a homogeneous, isotropic, elastic, multi-layered system, the calculated stress-strain distribution in the structure under the simple action of a passing wheel is very complex. For a realistic modeling of pavements, each layer and material constituent should be modeled with the appropriate material constitutive law, and the interaction of each layer should be analyzed as a pavement system. The finite element method is the most versatile approach in the mechanistic pavement design approaches. The FEM is an integral part of the PANDA approach; the modeling space (2D or 3D), wheel loading configuration, constitutive law, etc., are of great importance [67,69,76,79]. Collop et al. [57] implemented Burgers analogical constitutive model into FEM to predict viscoplastic deformation. A comprehensive FEM study by [80] investigated the effect of loading scenario in 2D and 3D modeling spaces and implemented the PANDA constitutive law. They investigated 11 different loading modes (given in Table 4) and compared three different material constitutive equations: elasto-viscoplastic, viscoelastic-viscoplastic, and viscoelastic-viscoplastic-viscodamage or PANDA. Based on wheel tracking test data, they concluded that the 2D model significantly overestimates permanent deformation (rutting) compared to the 3D moving load case (which is the most realistic case). Moreover, they found that the viscoelastic-viscoplastic-viscodamage model gives higher rutting damage predictions compared to elasto-viscoplastic and viscoelastic-viscoplastic models. The viscoelastic-viscoplastic and elasto-viscoplastic models give close predictions, where the form gives a still higher rutting prediction. On the contrary, another study using the viscoelastic-viscoplastic-viscodamage constitutive model and 3D FEM modeling of tire-pavement interaction showed that the pulsatile and equivalent loading assumptions overestimated rutting compared to the realistic moving load [90]. The FEM simulation has been applied for different viscoelastic, viscoplastic, and crack modeling. The linear viscoelastic simulation has been modeled using the analogical models (e.g., Maxwell model), the Prony series, and time-temperature superposition principles in commercial software such as ABAQUS.

Permanent Deformation and Fatigue Interaction Damage
Permanent deformation and fatigue cracking are the two dominant pavement damage mechanisms. Extensive research has been conducted to characterize the two damages independently. Advanced mechanistic models were proposed [64,83,84,91]. Fatigue damage occurred due to the formation and propagation of cracks due to repetitive loads. The bottom-up cracking of the asphalt concrete layer was the traditional fatigue mechanism. However, experimental and field observations showed that top-down cracking due to tire compression is also the cause of fatigue damage [92,93]. The crack initiation in compression occurs when the viscoplastic strain hardening reaches saturation at the flow number. The extra energy at maximum saturation (hardening) is consumed to initiate microcracks and increase phase angle. When asphalt with pre-existing cracks is subjected to a compressive load, wing cracks develop and propagate parallel to the load direction [94,95], contrary to tensile loading (cracks grow perpendicular to the stress direction). When pavement temperature is considered, fatigue cracking is critical at low temperatures, and permanent deformation is a high-temperature damage. At elevated temperatures, the critical energy threshold (the mixture relaxes faster) becomes more resistant to micro-cracks and needs more energy to initiate cracks [96]. On the other hand, deformation via plastic flow (aggregate re-orientation) is dominant at high temperatures. When asphalt concrete is subjected to repetitive loading, energy is dissipated by viscous flow and/or plastic flow, leading to fatigue cracking and/or permanent deformation [97], and some part of the energy is transferred into heat [50]. The energy dissipation caused material ductility exhaustion, hardening, and viscoplastic flow. The classic energy balance principle states that the decreasing rate of potential energy (stored and recoverable) during crack initiation is equal to the dissipated energy rate due to plastic/viscoplastic deformation and crack opening, and different failure criteria were proposed [98,99]. The dissipated energy (W) in a cyclic load is the area under the stress-strain hysteresis and expressed as follows.
where τ is the integration variable, σ is the stress function, and ε is the strain function. One hysteresis cycle consists of the total strain energy, which is the sum of damage-causing dissipated energy (the hysteresis area) and the recovery strain energy (stored energy). The classic approach considered rutting and fatigue separately (cf. Figure 15a) [100]. However, both fatigue and permanent deformation damages are caused by the same load in pavements. The difference is which damage mode dominates at different temperatures. Therefore, a realistic damage modeling should consider the simultaneous damage evolution of fatigue and permanent deformation. An interaction domain of rutting and fatigue can be shown in Figure 15b [101]. For an interaction damage condition, the dissipated energy (in Equation (26)) can be decomposed further into permanent deformation and fatigue damagecausing energies. As discussed above, the existing models for permanent deformation or fatigue damages are based on the idealized hysteresis loops of pure creep and pure fatigue. Only very few attempts can be found; for example, Luo et al. [99] investigated the permanent deformation behavior of pre-cracked asphalt concrete.
In a cyclic creep-recovery test, material hardening grows, and the hysteresis loops shift horizontally due to the accumulation of irrecoverable viscoplastic strain (open loops). On the other hand, in the idealized fatigue test (either stress-or strain-controlled), cyclic load is applied that causes stiffness reduction and a phase angle increment with a negligible viscoplastic strain (hysteresis loops do not shift horizontally or loops are closed). In simultaneous fatigue-permanent deformation damage, the material hardening contributes to initiate fatigue cracking at intermediate temperatures, and the fatigue cracks accelerate the microstructure change and contribute to viscoplastic deformation. The viscous flow creates a plastic zone (the crack initiation point), and fatigue cracking can evolve without plastic flow (aggregate movement and re-orientation). That means the simultaneous occurrence of the two damages cannot be ignored as they are interdependent. Thus, the hysteresis loops can be modified by superposing the two loops (Figures 4a and 16a) to account for the interaction of the two damages, as shown in Figure 16b. Similar fatiguecreep interaction loops can be found elsewhere for steel [102,103]. A detailed discussion of the interaction of the two damages is presented in next sections.
The classic approach considered rutting and fatigue separately (cf. Figure 15a) [100]. However, both fatigue and permanent deformation damages are caused by the same load in pavements. The difference is which damage mode dominates at different temperatures. Therefore, a realistic damage modeling should consider the simultaneous damage evolution of fatigue and permanent deformation. An interaction domain of rutting and fatigue can be shown in Figure 15b [101]. For an interaction damage condition, the dissipated energy (in Equation (26)) can be decomposed further into permanent deformation and fatigue damage-causing energies. As discussed above, the existing models for permanent deformation or fatigue damages are based on the idealized hysteresis loops of pure creep and pure fatigue. Only very few attempts can be found; for example, Luo et al. [99] investigated the permanent deformation behavior of pre-cracked asphalt concrete. In a cyclic creep-recovery test, material hardening grows, and the hysteresis loops shift horizontally due to the accumulation of irrecoverable viscoplastic strain (open loops). On the other hand, in the idealized fatigue test (either stress-or strain-controlled), cyclic load is applied that causes stiffness reduction and a phase angle increment with a negligible viscoplastic strain (hysteresis loops do not shift horizontally or loops are closed). In simultaneous fatigue-permanent deformation damage, the material hardening contributes to initiate fatigue cracking at intermediate temperatures, and the fatigue cracks accelerate the microstructure change and contribute to viscoplastic deformation. The viscous flow creates a plastic zone (the crack initiation point), and fatigue cracking can evolve without plastic flow (aggregate movement and re-orientation). That means the simultaneous occurrence of the two damages cannot be ignored as they are interdependent. Thus, the hysteresis loops can be modified by superposing the two loops (Figures 4a and 16a) to account for the interaction of the two damages, as shown in Figure 16b. Similar fatiguecreep interaction loops can be found elsewhere for steel [102,103]. A detailed discussion of the interaction of the two damages is presented in next sections. One of the limitations for interaction damage modeling is the lack of integrated testing protocol for permanent deformation and fatigue. Some attempts were made, such as Indirect Tension (IDT) testing for deformation and fracture and applying a haversine loading waveform in fatigue tests [104,105]. The haversine load is used with the assumption that the pulse can be decomposed into pure creep and pure sinusoidal components. However, the creep-recovery behavior cannot be captured with such approaches. In another study, Zhang et al. [106] conducted a destructive dynamic modulus test to simultaneously characterize permanent deformation and fracture properties using 16 different asphalt mixtures. In their work, they successfully modeled viscoelastic, viscoplastic, and viscofracture properties using a new destructive dynamic modulus test. Another reason for the independent treatment of fatigue and rutting damages on asphalt concrete damage is the perception that permanent deformation (creep) is expected in the early life of pavements, while fatigue is high-cycle damage after asphalt concrete accumulates sufficient hardening and age. However, this precedence is dependent on several factors such as temperature, loading range (loading time, rest period, deviatoric stress), material stiffness, number of cycles, and other climatic factors. For example, fatigue cracking can evolve before creep damage in stiff and thick pavements. The top-down cracking due to tire compression can develop without permanent deformation. Some field and laboratory observations have also shown that fatigue cracking accompanied rutting [92]. Therefore, both the creep-fatigue and fatigue-creep damage sequences can occur in asphalt concrete pavements [99].
The main objective in the design and service life estimation of asphalt concrete is to estimate the number of cycles (N ) to initiate fatigue cracking and/or the critical strain (ε - One of the limitations for interaction damage modeling is the lack of integrated testing protocol for permanent deformation and fatigue. Some attempts were made, such as Indirect Tension (IDT) testing for deformation and fracture and applying a haversine loading waveform in fatigue tests [104,105]. The haversine load is used with the assumption that the pulse can be decomposed into pure creep and pure sinusoidal components. However, the creep-recovery behavior cannot be captured with such approaches. In another study, Zhang et al. [106] conducted a destructive dynamic modulus test to simultaneously characterize permanent deformation and fracture properties using 16 different asphalt mixtures. In their work, they successfully modeled viscoelastic, viscoplastic, and viscofracture properties using a new destructive dynamic modulus test. Another reason for the independent treatment of fatigue and rutting damages on asphalt concrete damage is the perception that permanent deformation (creep) is expected in the early life of pavements, while fatigue is high-cycle damage after asphalt concrete accumulates sufficient hardening and age. However, this precedence is dependent on several factors such as temperature, loading range (loading time, rest period, deviatoric stress), material stiffness, number of cycles, and other climatic factors. For example, fatigue cracking can evolve before creep damage in stiff and thick pavements. The top-down cracking due to tire compression can develop without permanent deformation. Some field and laboratory observations have also shown that fatigue cracking accompanied rutting [92]. Therefore, both the creep-fatigue and fatigue-creep damage sequences can occur in asphalt concrete pavements [99].
The main objective in the design and service life estimation of asphalt concrete is to estimate the number of cycles (N f ) to initiate fatigue cracking and/or the critical strain (ε c -rupture strain) that causes the flow or rutting of asphalt concrete. The simplest way to evaluate the interaction damage of creep and fatigue is separately calculating the creep and fatigue damages and adding them together.
By the life fraction rule, the sum of fatigue damage (φ f ) and creep damage (φ c ) equals a certain damage density value, φ. The nonlinear summative decomposition of damage can be expressed as follows for a more general representation.
The classic continuum damage models were used to define the damage rates for creep [63,65] and fatigue damage [107,108]. Lemaitre et al. [109] proposed sequential damage interaction, and the total damage accumulation during one cycle of creep followed by fatigue and fatigue followed by a creep sequence can be expressed as follows (respectively).
For time-dependent materials such as asphalt concrete, damage densities φ c_f and φ f_c will not be the same for the same number of cycles in both sequences. This is because of the different modes of damage formation in creep and fatigue and the possible interactive damage one on another. A parameter can be defined for "interactive" damage. Skelton et al. [101] presented analytical expressions to allow creep to be modified by fatigue and fatigue to be modified by creep. The combined equation for the "creep-fatigue" and "fatigue-creep" interactive damages can take the following form.
The interaction coefficients I cf (creep on fatigue) and I fc (fatigue on creep) can take any value between zero and unity. The damage density for fatigue and creep can be determined by rearranging terms and solving a quadratic solution. In the creep-fatigue sequence, it is assumed that pure creep is followed by fatigue damage. The fatigue damage part is modified to account for the pre-existing creep damage. Therefore, the creep-fatigue and fatigue-creep interactive damages can be expressed as follows, respectively.
Re-arranging and solving for φ c and φ f gives the following expressions.
The creep-fatigue and fatigue-creep interaction coefficients I cf and I fc are non-zero and take any value between −1 and 1. Figure 17 shows the relationship between creep and fatigue damage densities at random values of creep-fatigue interaction coefficients.
Re-arranging and solving for ϕ and ϕ gives the following expressions.
The creep-fatigue and fatigue-creep interaction coefficients I and I are non-zero and take any value between −1 and 1. Figure 17 shows the relationship between creep and fatigue damage densities at random values of creep-fatigue interaction coefficients.

Summary
The practice is always trailing the theory in pavement design, and the case of the permanent deformation prediction is no different. Many highway agencies and road administrations are practicing an empirical method of permanent deformation characterizations (for example, the Marshall test). Some progresses are being made to transition from empirical to mechanistic-empirical design methods. The primary motivation of this literature study was to explore the advancement of asphalt concrete permanent deformation characterization, constitutive modeling, and application. The current state of research is focused on the formulation of a mechanistic method that applies the fundamental theories of mechanics and materials to predict permanent deformation damage. In the last decade, promising advancements have been made in the development of comprehensive and coupled permanent deformation modeling. Pavement analysis using the nonlinear damage approach (also known as PANDA) constitutive model is one of the notable progresses that is gaining wide acceptance (as the next generation mechanistic pavement design approach). Another promising area of progress is the microstructural approach aided with 3D digital image analysis equipment and finite element modeling. This digital technology-based permanent deformation modeling is also the future prospect for accurate permanent deformation modeling and prediction of different asphalt mixtures. The mechanistic methods offer unlimited potential to expand the modeling parameters, and different damages and phenomena such as aging, healing, moisture damage, pre-crack, etc., can be coupled to unify damage prediction. The drawbacks of mechanistic models are sophistication, requiring extensive testing, and the calibration of several modeling variables. For example, the viscoelastic-viscoplastic-viscodamage model requires more than 21 model variables to be optimized using at least two different experiments and several test repetitions at different temperatures stresses and strains. The improvement of the latest models from the classic viscoplastic strain hardening model is the consideration of cyclic hardening and relaxation mechanisms and the viscodamage of asphalt concrete. Although the mechanistic method is theoretically appealing, the calibration cost and rigorous equations can be considered the limitations. The micromechanics approach considers the evolution of permanent deformation related to changes in the microstructure of asphalt concrete constituents (aggregates and mastic). Thus, the micromechanics method is regarded as the most realistic way of modeling heterogeneous materials such as asphalt concrete. Based on the extensive literature study, the permanent deformation prediction and modeling approaches can be categorized into four aspects, as shown in Table 5.

Conclusions
A review study was conducted to explore the state of the art on permanent deformation prediction from the 1960s' pure empirical to latest mechanistic (from 2011) methods. The review study revealed that the latest constitutive models integrate and couple different theories, i.e., continuum mechanics (1958), nonlinear viscoelastic (1969), viscoplasticity (1971), and viscodamage (2011), along with the crucial time-temperature superposition principle. Such coupling techniques offered advantages of integrating different asphalt concrete damages and opened the possibility of a unified asphalt damage model in the future. The PANDA model is one of the most comprehensive permanent deformation modeling approaches available in the literature (the next generation to the mechanistic-empirical method). The calibration and/or validation tests are reliant upon the conventional creep or creep-recovery tests in either confined or unconfined modes. The computation and experiment cost of mechanistic methods are the limitations. The practical application of the mechanistic models is very limited at the moment. Moreover, the fatigue-permanent deformation (rutting) interaction is often ignored in the existing (studied) literature. It is inferred that both damages can evolve simultaneously as the same load caused both damages. The mechanistic approach has the potential to couple the two predominant pavement damages. From the extensive study, it can be synthesized that a unified permanent deformation damage model can be developed by integrating continuum damage and microstructure approaches and coupling fatigue, moisture, healing, aging, and other physical and chemical phenomena in asphalt concrete.

Further Research
The mechanistic method is "universally" applicable to predict damage regardless of climatic conditions, stress state, or material type. This characteristic presents wide, open research questions, for example, (1) developing a unified pavement damage performance prediction model, (2) a coupled model for fatigue and rutting damages using continuum mechanics and viscoelastic and viscoplastic theories, (3) considering (coupling) different strains such as shear and axial strains in permanent deformation prediction models, and (4) developing simplified (unified) asphalt concrete test methods to characterized different damages simultaneously (fatigue, rutting, moisture, etc.). The authors of this paper are conducting research on the simultaneous creep-fatigue damage evolution in a sequential manner.