Continuum-based micromechanical models for asphalt materials: Current practices & beyond

The mechanical properties of asphalt mixture are always required for the evaluation of the durability of pavements. In order to obtain these properties without conducting expensive laboratory tests and using calibrated empirical models, research studies have been carried out to develop micromechanics-based models. Continuum-based micromechanical models (CBMM), which are developed based on continuum mechanics, have increasingly been utilized to estimate the mechanical properties of asphalt materials based on the fundamental properties of individual constituents. These analytical models are expected to provide reliable predictions without the need for extensive computational facilities. Although the utilization of CBMM has been presented by several past studies, most of the studies do not provide a concise and critical review of these models. Therefore, in this paper, a complete review of CBMM was presented. Commonly used CBMMwere introduced in detail and their advantages and disadvantages were discussed and compared. Comprehensive summaries and critical discussions about their current utilization and limitations for predicting the mechanical properties of asphalt materials were given. Further modifications and new development for addressing the limitations were extensively described and discussed. In the end, research challenges were highlighted and future recommendations from different perspectives were proposed. 2020 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Asphalt mixture is the most commonly used material in pavement construction. Pavement engineers/researchers generally measure the mechanical properties of asphalt mixture such as dynamic modulus, creep compliance, etc. to evaluate the durability of pavements in the field. Therefore, in commonly used pavement design tools such as Mechanistic-Empirical Pavement Design Guide (MEPDG), these properties are often required as material inputs. In order to obtain the mechanical properties of an asphalt mixture, laboratory tests are usually conducted. On one side, these tests are reliable; however, on the other side, they require a series of expensive sampling and laboratory testing facilities. Moreover, some tests require a long time to perform.
As an alternative to laboratory tests, researchers have developed different empirical models such as the Heukelom and Klomp model [1], the Bonnaure model [2], Witzak's model [3,4], etc. These models have been effectively utilized in obtaining the mechanical properties of asphalt mixtures on the basis of the properties of aggregates and binder. The drawbacks of such empirical models are that: a) they are developed by statistical regression analyses which do not give much insight into the fundamental mechanisms behind mixtures' behaviors; b) they are only applicable on the mix-ture which they are calibrated for, and when they are used for different mixtures, the accuracies of the models decrease.
To deal with the above-mentioned limitations of empirical models, research studies have been carried out to develop socalled ''micromechanics-based models". These models evaluate the localized stress and strain fields for a given macroscopic loading condition. The obtained stress/strain fields are further utilized to obtain a composite's effective properties with the fundamental properties of individual constituents.
Numerical micromechanical models, i.e. finite element models (FEM) and discrete element models (DEM), have been successfully utilized by many researchers [5][6][7][8][9][10][11]. These models can handle complex compositions of a mix by capturing its realistic microstructure characteristics using the Computed Tomography (CT) scan technique [12]. However, the CT scan technique still requires dedicated equipment and post-processing facilities, which are not readily available to pavement engineers. Recent studies [13,14] have shown that a random assembly of aggregates can also be virtually generated using an image-aided random generation method based on a pre-built aggregate image database. This process still requires post-processing tools for FEM/DEM-meshes creation and, moreover, simulations of a mix's properties require large-scale computational facilities. Since such computational tools and facilities are not typically available to the field engineers, such models cannot be widely used in practice. Therefore, analytical micromechanical models are seen as a good alternative which is expected to provide reliable estimations of mechanical properties of a composite without the need for extensive computational facilities.
Analytical micromechanical models developed based on continuum mechanics, continuum-based micromechanical models (CBMM), have increasingly been used to predict the mechanical properties of asphalt materials. In such models, the detailed information of individual constituents is not required to be described. On the contrary, the constituents having the same (very similar) mechanical properties are treated as one phase; and a composition of various phases represents a composite. For a given macroscopic loading condition, each phase's average stress and strain are evaluated and further utilized to obtain the effective properties of the composite on the basis of the volumetric, mechanical and/or geometrical properties of individual phases.

Application of CBMM for asphalt mixture
Using CBMM, the effective mechanical properties of an asphalt mixture are determined on the basis of individual phases' properties such as aggregates' stiffness, binder's stiffness, volume fractions, etc. The obtained effective properties are further used to calculate the local stress/strain fields (at the component level) for predicting the mixture's performance under a given loading condition [15][16][17], see Fig. 1.
From the pavement engineering perspective, CBMM can be seen as an alternative to handle the following research needs: a. evaluation of the stresses and strains responses of asphalt mixture at critical locations in the design of pavement. This can be of particular importance for unconventional mixtures. b. prediction of the propensity of a given mixture to different distresses. For example, the predicted local stress/strain fields can be utilized to analyze raveling of porous asphalt mixture by evaluating the fatigue life of the mastic and the mastic/aggregates interface, Fig. 1. c. optimization of source materials composition to achieve desirable mechanical properties of a mixture.
Several past studies have presented the utilization of CBMM. However, most of the studies do not present a complete overview of these models in detail. Moreover, a comprehensive description of their current utilization with/without modification is also not available to the researchers and practitioners. Considering the significant possibility of these models' utilization for asphalt materials, this paper seeks to provide a concise but critical review of these models.

Scope of the research study
The scope of this study includes the following: future recommendations for further improving the performance of CBMM.

Introduction of CBMM
2.1. Definition of the effective stiffness of a composite 2.1.1. Stiffness tensor In linear elasticity, the constitutive law of a material is generally given as: where r is the second-order stress tenser; e is the second-order stain tenser; C is the fourth-order stiffness tensor; the symbol ':' means the double dot product between two tensors.
For isotropic elastic materials, five material constants, i.e. Young's modulus E, shear modulus G, bulk modulus K, Lame constant k, and Poisson's ratio m, are commonly used. However, only two of them are required to describe the characteristic of C completely [18]. For example, by using the values of K and G, C can be represented as: In Equation (2), I v and I d denote the volumetric part and the deviatoric part of a four-order tensor, respectively; and they are defined as where d is the Kronecker's delta. From these definitions, the following relations between I v and I d can be derived: where I is the unit fourth-order tensor, which is written as On the basis the relations in Equation (4), it can be found that if two fourth-order tensors B 1 and B 2 are represented as the double dot product of B 1 and B 2 can be directly given by 2.1.2. Effective stiffness of a composite According to Equation (1), in CBMM, the effective stiffness of an N-phase composite C eff is defined using the composite's average stress <r> c and average strain <e> c : The values of <r> c and <e> c for a given volume V of a representative volume element (RVE) can be described by Equation (10).
At the component level, the average stress <r> r and strain <e> r of phase r over the volume of this phase (V r ) are given as Equation (11).
Substituting Equation (11) into Equation (10), <r> c and <e> c are rearranged as where / r is the volume fraction of phase r in the RVE, which is defined as For each phase, it is also known that <r> r and <e> r satisfy the constitutive law: where C r is the stiffness tensor of phase r. By substituting Equation (9) and (15) into Equation (12), the value of C eff can be related to the values of C r : Equation (17) shows a typical relationship between <e> r and <e> c which is used in CBMM [19]: where A r is called the strain localization tensor of phase r. By combining Equations (16)and (17), the value of C eff is further expressed as: Equation (19) can be easily derived once Equation (17) is substituted into Equation (13).
Equations (18) and (19) show that the values of / r , C r and A r must be determined or estimated to obtain the value of C eff . In general, the constituents of a composite are known, which means the values of / r and C r can be either determined in the laboratory or readily available in the literature. It is highlighted here that although the value of A r , to a certain extent, can also be measured by sophisticated technologies such as digital image processing, smart sensors, etc., they are not frequently available to pavement engineers and researchers. The central idea behind CBMM is primarily to calculate the value of A r .

Approaches for predicting a composite's effective stiffness
Mainly two different approaches are used to estimate a composite's effective stiffness. A brief description of these approaches is given in the following paragraphs.

Geometry based approach
Geometry based approach is based on a geometric model [20] in which the relative arrangements of individual phases for a given composite are predefined, and in general, a closed-form solution for the effective stiffness can be obtained. This approach can be further subdivided into following subcategories: a) models developed from an arrangement of individual phases in parallel or series or a combination of them (see Fig. 2a); for examples, the Voigt model [21]; the Reuss model [22] and the Hirsch model [23], and b) models developed by describing a composite as a matrix with different inclusions embedded into it (see Fig. 2b); for examples, the Dilute model, the Self-consistent (SC) model, the Generalized selfconsistent (GSC) model, the (N + 1)-phase model, the Mori-Tanaka (MT) model, the Differential model, and the Composite sphere model, etc.
The solutions for effective moduli provided by phase arrangements (i.e. subcategory a) are easy to implement; however, they most often require calibration factors. Whereas, Eshelby-based micromechanical models (i.e. subcategory b) are developed on the basis of the Eshelby's solution [19] for the inhomogeneity problem where an ellipsoid inclusion is embedded into an infinite matrix. By using the Eshelby's solution, the relationship between the effective stiffness of a composite with the properties of each phase can be obtained.

Bounds based approach
The other approach to estimate the effective stiffness of a composite is to determine the lower and upper bounds of the possible range of the stiffness [24][25][26][27][28][29][30][31]. In the derivation of these bounds, no specific geometry of the composite is defined, and only limited available information is utilized. For example, by using only the volume fraction and mechanical properties of each phase, Paul's bounds [24] and Hashin-Shtrikman (H-S) bounds [25,27], are derived. The H-S bounds are considered as the most restrictive bounds which can be given in terms of phases' moduli and volume fractions [25,27]. The derivation of more restrictive bounds requires additional information about the structure of the composite [30,[32][33][34].
It is noted that in addition to estimating the effective stiffness, the bounds based approach is often used to validate CBMM. It is considered that a model is possible to perform well only when its predictions lie between the upper and lower bounds.

Commonly used CBMM
This section presents complete descriptions of the commonly used CBMM developed using the geometry-based approach, i.e. the Voigt model, the Reuss model, and the Eshelby-based micromechanical models, and the bounds based approach, i.e. Paul's bounds and the H-S bounds.
By combining Equation (16) and Equation (20), the relationship between C eff and C r can be obtained as When the composite is considered to be isotropic, Equation (21) is expressed using the effective bulk modulus K eff and the effective shear modulus G eff as where K r and G r denote the bulk modulus and the shear modulus of phase r, respectively. On the contrary, in the Reuss model, all the phases are assumed to be in a series arrangement, Fig. 4. The stresses of individual phases are the same, Equation (23). The relationship between C eff and C r can be obtained by substituting Equation (15) and Equation (23) into Equation (13), see Equation (24). By considering the composite as isotropic, the values of K eff and G eff can be written as Equation (25).
The Voigt model and the Reuss model use simple assumptions to relate the stress/strain filed of each phase to the stress/strain field of the composite. They can be easily implemented; however, the calculated results of C eff are generally inaccurate [35,36]. To solve this issue, researchers have proposed more complicated arrangements of individual phases. For example, for predicting concrete's modulus from the properties of cement and aggregates, Hirsh [23] proposed a model where the cement phase and the aggregate phase are arranged in a combination of parallel and series, as shown in Fig. 2a. The total volume of each phase is separated into the series portion and the parallel portion. Using the same method, Christensen et al. [37] proposed an arrangement to predict the modulus of asphalt mixture. The predictions of these models were observed to be in good agreement with the experimental values. However, the calibration of the proportion of the series and the parallel part of each phase is always required for accurate predictions [16].
where C 1 and C 2 denote the stiffness tensors of the matrix and the inclusion, respectively; and S 1 is known as the Eshelby's fourthorder tensor.
The value of S 1 is a function of the matrix's mechanical properties and the inclusion's mechanical and geometrical properties [38][39][40][41]. When the matrix and the inclusion are isotropic materials and the inclusion is a sphere, S 1 is calculated as with where K 1 denotes the matrix's bulk modulus, and G 1 denotes the matrix's shear modulus.  (27) with C r , see Equation (29). The value of A 1 is further obtained by using Equation (19).
When the inclusions are considered as spheres and all the phases are isotropic materials, the values of K eff and G eff are given as Since the values of A r are obtained directly from Eshelby's solution in which the matrix phase is considered as infinite, the Dilute model is only suitable for composites in which the inclusions' concentrations are low enough to neglect the interaction between stress/strain fields disturbed by different inclusion particles (known as ''inter-particle interaction").
With the aim of improving the predictions' accuracy, researchers have proposed the self-consistent scheme which considers the inter-particle interaction. In this scheme, either the inclusion itself (the SC model) or the inclusion together with its surrounding/s (the GSC/(N + 1)-phase model) is considered to be included into an infinite medium with identical characteristics as the composite itself. In the following section, the models using the concept of the self-consistent scheme are presented.
In the SC model, the value of A r for each inclusion phase is calculated by replacing the mechanical properties of the matrix in Equation (27) with the unknown mechanical properties of the composite (C eff ), see Equation (32). The value of A 1 can still be obtained from Equation (19).
A r ¼ ½I þ S eff : ðC eff Þ À1 : ðC r À C eff Þ À1 ; r ¼ 2; 3; :::; N ð32Þ K eff and G eff for isotropic and spherical inclusions can be calculated using Equations (33) and (34). The expressions of K eff and G eff are not explicit, and thus numerical techniques should be used to solve these equations.
It is noted that in the SC model, each phase is treated equally [46]. Therefore, the predictions do not depend on the selection of the matrix and the inclusion, while they are controlled by the phase which has the highest volume fraction. For example, for asphalt mixture, which typically has a volume fraction of aggregates of more than 50%, the predicted effective stiffness is more similar to the properties of aggregates which are stiff and frequency/temperature independent [15].  [47], see Fig. 7. It is noted here that although it seems that the model contains three phases, in reality, it is only suitable for a two-phase composite as the infinite medium is the composite itself.
Regarding an isotropic and spherical inclusion, the solutions for the values of K eff and G eff of a two-phase composite are given as Equation (35) and (36), respectively.
Coefficients of A, B and C are calculated using the following equations: The GSC model has been found to give more physically reasonable predictions of the inter-particle interaction for a two-phase composite in comparison to the MT model and the Differential model [48]. However, it is highlighted here that the GSC model does not provide an analytical solution for the effective shear modulus of a multi-phase composite. Therefore, a sequence step method is generally used where only one inclusion phase is considered in each step [49,50]. For example, to predict the effective stiffness of asphalt mixture, the effective stiffness of the mix of asphalt binder and aggregates (or air voids) is firstly predicted. After that, by considering the obtained mix as the matrix, the air voids (or aggregates) are considered as inclusions to calculate the mixture's effective stiffness. The limitation of this method is that the predictions depend on the order of adding different phases, especially when the volume fractions of the inclusions are similar to each other, such as a porous asphalt mixture with a high volume fraction of air voids [15].
2.3.1.6. The (N + 1)-phase model. By adding more layers surrounding the central spherical inclusion, Herve and Zaoui [51] extended the GSC model to the (N + 1)-phase model, Fig. 8. In the figure, the central sphere with a radius of R N is composed of phase N; and the shell formed by spheres with radii of R r+1 and R r is composed of phase r.
By considering all phases as homogeneous and isotropic, the value of K eff can be computed from Equation (37); and the value of G eff can be obtained by solving a quadratic formula, Equation (38). It is highlighted here that when N is equal to two, these equations are identical to those in the GSC model (Equation (35) and (36)).
Coefficients of A, B and C can be obtained from the following equations: where m r is the Poisson's ratio of phase r, and R r is the radius of each phase, which is calculated from the volume fraction of each phase (/ r ) using Equation (39).
The (N + 1)-phase model is generally utilized when an additional phase is required to be modeled between the inclusion and the matrix. For instance, this phase can be an actual coating material or result from physiochemical interactions between different phases [51].
2.3.1.7. MT model. In the pioneering study of Mori and Tanaka [52], an approach to obtain the matrix's average stress was described. In later studies, researchers [53][54][55][56][57][58] further utilized the concept of average stress together with Eshelby's solution (as discussed before) to obtain the effective properties of a given composite.
The MT model assumes that inclusions are included into a medium which processes the same properties as the matrix, Fig. 9, and it is assumed that the value of <e> r for each inclusion phase is calculated from Equation (40).
where T r is identical to T by replacing C 2 with C r in Equation (27). By substituting Equation (40) into Equation (13), <e> 1 is related to <e> c via Equation (41). From Equation (41), the value of A 1 is directly obtained as Equation (42); whereas the value of A r for each inclusion phase is obtained by combining Equations (40) and (41), see Equation (43).
The values of K eff and G eff for isotropic and spherical inclusions can be calculated from Equation (44) and (45), respectively. It is highlighted that when N is equal to two, the solution for K eff in Equation (44) is the same as that in Equation (35). Fig. 9. Illustration for the MT model.

R N R r+1
Phase N=inclusion

Phase r+1
Phase r R r R 1

Phase 1
Infinite medium Comparing to the self-consistent scheme, the calculation of K eff and G eff via the MT model is easier to be implemented since no implicit equations and quadratic formula are required to be solved. However, it was found by researchers [48] that the inter-particle interaction approximated by the MT model is less accurate than the self-consistent scheme, especially under concentrated conditions.
2.3.1.8. The Differential model. The Differential model is another method to deal with the inter-particle interaction in a highly concentrated composite. The idea behind this model is to develop a process where the interactions can be neglected. This model was initially developed for fluid suspensions [59]. McLaughlin [60] further showed that this model can also be utilized in the field of continuum micromechanics. Norris [61,62] provided a more generalized version of the model to handle multiphase composites. In the following part, the derivation of the Differential model is explained by taking a two-phase composite as an example, Fig. 10.
In the first step, a small number of inclusion particles V 2 1 are added into the matrix with a volume of V 1 to obtain a composite which is called ''Effective medium 1". Based on Equation (18) and (19), the effective stiffness of this medium C eff 1 can be calculated as: Since the embedded particles are in dilute condition, the interparticle interaction can be neglected, and the value of A 2 0 can be calculated via the Dilute model, Equation (47).
The process is repeated by treating ''Effective medium 1" as a new matrix and adding another small volume of inclusion particles to obtain ''Effective medium 2" with a stiffness tensor of C eff 2 . This iteration is continued until the total volume fraction of the inclusion is the same as that in the composite. In the r th step, if the added inclusion's volume is V 2 r , the value of C eff of the obtained ''Effective medium r" is given by: where A 2 r-1 is calculated from the properties of the effective medium obtained in step r-1: The total volume fraction of the inclusion in the r th step, / r 2 , is written as and in the previous step, / rÀ1 2 is given by From the rearrangement of Equation (50) and (51), the following equation is obtained: where D/ r 2 = / r 2 -/ rÀ1 2 . By substituting Equation (52) into Equation (48), Equation (53) can be obtained: From Equation (53), an approximated value of C eff can be obtained by an iterative process with a small value of D/ r 2 . Considering the value of D/ r 2 to be infinitesimally small, Equation (53) in a differential form can be written as Equation (54). where Equation (54) can be solved by separating the variables and then integrating them on both sides: If C is represented in terms of K and G, Equation (56) can be rewritten as where A 2 v and A 2 d are the volumetric and deviatoric parts of the fourorder tensor A 2 . Assuming both phases as isotropic materials and the inclusion particles as spheres, Equation (57) and (58) can be further expressed as: Effective medium 2 Fig. 10. Illustration of the Differential model.
Equation (54) for a three-phase composite can be expressed as [61]: where / 2 ðtÞ and / 3 ðtÞ are the volume fractions of phase 2 and phase 3 at every step ''t", respectively, and / t ð Þ ¼ / 2 t ð Þ þ / 3 ðtÞ. In order to obtain a unique value of C eff from Equation (61), a path must be defined from zero to the final values of / 2 and / 3 in the (/ 2 , / 3 ) plane considering that the value of C eff is pathdependent [61,62]. By maintaining the volume fractions of inclusions to be constant in each step (see Equation (62)), Equation (61) can be rearranged as Equation (63). It is highlighted here that in Equation (63), ''t" is eliminated and thus the value of C eff can be obtained.
where the superscript ''c" indicates the final composite.
Even though a specific geometry is defined in the Composite spheres model, only the exact solution for the value of K eff , which is the same as the solution in the MT model (Equation (44) G eff U ) and the lower bounds (K eff L and G eff L ) are given as the predicted moduli using the Voigt model (Equation (65)) and the Reuss model (Equation (66)), respectively.
Although Paul's bounds are easy to be obtained, these bounds are generally not accurate enough to provide good estimates of the effective moduli [35,36]. For example, for a two-phase composite, researchers [35,36] concluded that when the ratio of these two phases' moduli is higher than 2:1, the distance between upper bound and lower bound becomes too wide for any practical utilization.

The H-S bounds.
The H-S bounds are developed by using variational principles [25]. The general solutions for the H-S bounds proposed by Walpole [27] are shown in Equations (67)- (70).
À1 where K max and K min represent the maximum and minimum values of the bulk modulus among all phases of the composite, respectively; and G max and G min denote the maximum and minimum values of the shear modulus, respectively. Although the H-S bounds can provide better estimates for the effective moduli of a composite in comparison to Paul's bounds, these bounds still fail to provide narrow band results when the differences of moduli between individual phases are significant [27]. In summary, various CBMM have been developed to predict a composite's effective stiffness. The advantages and disadvantages of each model are summarized in Table 1. Depending upon the characteristics of the modelled composite, different models can be selected.

General solution procedure for viscoelastic composites
All the CBMM described above are originally developed for elastic composites; however, asphalt materials are mostly treated as viscoelastic composites. This means that these models may not be directly applicable to asphalt materials. Therefore, this section describes a typical procedure that can be adapted to utilize the above models for viscoelastic composites.
According to the research work of Hashin [64,65], micromechanical models can be utilized for viscoelastic materials via the elastic-viscoelastic correspondence principle [66]. Since the viscoelastic properties of a material are generally measured in the frequency domain and in the time domain, the general solution procedure of utilizing the models will be described in both domains. In the following part, a two-phase composite is taken as an example, and it is assumed that the matrix (phase 1) is a viscoelastic material while the inclusion (phase 2) is an elastic material. The equations of the Dilute model (Equations (30) and (31) In the frequency domain, the complex bulk modulus K 1 *(x) and the complex shear modulus G 1 *(x) of phase 1 can directly replace K 1 and G 1 , respectively. Since phase 2 is an elastic material, the moduli of this phase (K 2 and G 2 ) are left unchanged, Equations (71) and (72).
where K eff * and G eff * are the complex bulk and shear moduli of the composite, respectively. The calculated values of K eff * and G eff * can be further represented in terms of dynamic moduli (the absolute values of complex moduli) and phase angle. It is noted here that in some research studies, dynamic moduli are taken as the corresponding elastic moduli to be used in the formulas of micromechanical models [67][68][69]. However, this method can only obtain the effective dynamic moduli but not the phase angle of the composite.
In the time domain, the relaxation moduli or creep compliances of the viscoelastic phases are required to be transformed into the Laplace-Carson (LC) domain. For example, if the shear relaxation modulus R 1 d (t) (or the shear creep compliance D 1 d (t)) and the bulk relaxation modulus R 1 v (t) (or the bulk creep compliance D 1 v (t)) of phase 1 are described by the Maxwell model, Equations (73) and (74), their transformed formats in the LC domain are given as Equations (75) and (76).
where M and g are the spring's modulus and the dashpot's viscosity in the Maxwell model, respectively.
where s is a variable in the LC domain.
The effective relaxation moduli of the composite in the LC domain, R eff v* and R eff d* , can be obtained by replacing K 1 and G 1 in Equations (30) and (31)   Calibrations are required.

Dilute model Direct implementation of Eshelby's solution.
Cannot consider the inter-particle interaction.

SC model
The predictions are governed by the phase which has the highest volume fraction rather than the selection of the matrix and the inclusion. Consider the inter-particle interaction Implicit equations require to be solved.
GSC model Give more physically reasonable predictions of the inter-particle interaction for a two-phase composite; Analytical solution for effective shear modulus of a multi-phase (N > 2) composite is unavailable.
Can model additional phases between the inclusion and the matrix; Consider the inter-particle interaction.
Complex quadratic formula is required to be solved.
MT model Easily to be implemented; Consider the inter-particle interaction.
The approximated inter-particle interaction is less accurate. Differential model The inter-particle interaction can be ignored in each step. Complex differential equations are required to be solved. Composite sphere model Provide a more accurate geometry of a composite.
Exact solution for shear modulus is unavailable.

Bound based approach
Paul's bounds Easily to be implemented. Not accurate enough to provide good estimates of the effective moduli.

H-S bounds
The most restrict bounds which can be given in terms of phases' moduli and volume fractions.
Fail to provide narrow band results when the differences of moduli between individual phases are significant. Until now, the basic theories of CBMM and the general solution procedure for the utilization of these models in viscoelasticity have been presented. In later sections, the utilization of these models for predicting the viscoelastic properties of asphalt materials is presented and summarized. It is noted that the following summaries mainly focus on the Eshelby-based micromechanical models because these models do not require calibrations, and they can provide exact solutions for the moduli of asphalt materials.

Different scales of asphalt materials and upscaling
In general, a whole asphalt mixture can be subdivided into four length scales, i.e. the asphalt binder scale, the mastic scale, the mortar scale and the mixture scale, see Fig. 12 [49]. Although asphalt binder is a heterogeneous material as well, most research studies focus on the micromechanical modelling of the other three scales considering the fact that it is generally easy to directly measure the properties of asphalt binder from laboratory tests. Therefore, in this paper, asphalt binder is considered as a homogeneous material and the discussions mainly focus on the application of CBMM for other three scales.
The definitions of mastic scale, mortar scale and mixture scale are given as follows.
Mastic scale. Asphalt mastic is the material that results from the combination of asphalt binder and filler particles. It is generally assumed that air voids do not exist in asphalt mastic when the filler's concentration is not quite high (i.e. approximately less than 50% [70]). Hence, asphalt mastic is generally considered as a two-phase composite with filler embedded into asphalt binder [68,69,[71][72][73]. Mortar scale. Asphalt mortar is composed of asphalt binder, filler and fine aggregates. Until now, there is no consistent agreement about the properties of mortar in terms of the gradation of fine aggregates, the content of asphalt binder and the air voids content. Researchers [49,[74][75][76][77][78][79] prepared mortar samples using different sizes of fine aggregates, i.e. particles smaller than 2.36 mm, 2 mm, 1.18 mm, etc. The air voids content can vary from lower than 1% to higher than 20% [79] according to different aggregate gradations and asphalt binder contents. Asphalt mixture scale. Asphalt mixture is created by mixing asphalt binder, filler and graded aggregate particles. Generally, there are some air voids in an asphalt mixture. According to different design requirements, the air voids content can be lower than 5% (i.e. dense asphalt mixture) or higher than 20% (i.e. porous asphalt mixture).
The process of predicting the properties of a higher scale material from those of a lower scale material is known as upscaling. Upscaling can be done in one step to obtain the properties of a higher scale material directly. For example, to predict a mixture's stiffness, one way is to upscale the stiffness of asphalt binder by adding all the volumes of aggregates and air voids simultaneously.
Upscaling can also be conducted in a series of steps, which is known as multiscale upscaling [80][81][82][83]. In each step of the multiscale upscaling, the predictions from a lower scale are considered as the inputs of the matrix phase of an upper scale. By using this technique, the prediction of a mixture's stiffness can also be carried out as follows. In the first step, asphalt binder is considered as the matrix phase to predict the mastic's properties with the addition of filler as the inclusion phase. Then, the properties of mortar are predicted by considering the mastic obtained in the previous step as the matrix phase and the sand as the inclusion phase. Lastly, on the basis of the predicted properties of mortar as the matrix phase and the properties of stone and air voids as inclusions, the properties of the mixture are obtained.
The upscaling of asphalt materials can be conducted by using all the above introduced CBMM. Depending on the volume concentration of the inclusions, the test temperature, the order of the upscaling scale, etc., the performance of these models varies. In the following section, the sensitivities of the models on those factors are summarized.

Sensitivity of the models on test temperatures & volume concentrations
The performance of different models at different volume concentrations and test temperatures can be summarized as follows.
At very low concentrations, most of the models can obtain good predictions [69,71,72]. This is because all these models are developed on the basis of Eshelby's solution which is suitable for a composite with a low concentration of inclusions.

Asphalt binder
Mastic Mortar Mixture Upscaling Low High At high concentrations and low temperatures, the accuracy of the predicted results varies from one model to the other [49,70,73,81]. The Dilute model and the MT model have been generally found to under-predict the moduli (or over predict the creep compliances) of an asphalt material [79,81]. On the contrary, the SC model, the GSC model and the Differential model have been found to be more suitable for high concentrations of asphalt materials [15,79,81]. At high concentrations and high temperatures, none of these models have been found to provide accurate predictions, and in general, the predicted moduli are much lower than the measured values [69][70][71].

Sensitivity of the models on the order of the upscaling scale
Researchers [79,84] have found that upscaling from the properties of a higher scale matrix is more accurate than upscaling from a lower scale matrix. In some cases, the modulus of an asphalt material at high concentrations and high temperatures can even be accurately predicted by considering a higher scale material as the matrix [79]. This can be explained by the following two facts: In comparison to a lower scale matrix, more aggregate particles are included in a higher scale matrix. Therefore, the concentration of inclusions decreases, and thus the accuracy of the predictions improves. When a higher scale matrix is used for upscaling, the inaccuracy in the prediction from the lower scale matrix to the higher scale matrix is avoided because the properties of the higher scale matrix are accurately measured from laboratory tests.

Sensitivity of the models on the multiscale modeling technique
The multiscale modeling technique has been found to improve the predictions' accuracy in comparison to the upscaling conducted in one step [80,85]. To some extent, this technique can increase the accuracy of the predictions because the volume fraction of inclusions in each step is relatively low. Multiscale modeling has been successfully applied in particular testing conditions, such as low temperatures, while at high temperatures, the improvement has not been found to be very significant [15].

Limitations of the commonly used CBMM
The above section shows that the commonly used CBMM fail to accurately predict the properties of highly concentrated asphalt materials at high temperatures. In order to explain the limitation of these models, it is necessary to understand the mechanisms behind the stiffening of an asphalt material due to the addition of inclusions (known as ''stiffening mechanisms"). In this section, the stiffening mechanisms of asphalt materials are discussed, on the basis of which, the limitations of the models are explained.

Stiffening mechanisms of asphalt materials
There are three generally accepted stiffening mechanisms for asphalt materials [71]: the volume-filling reinforcement, the physiochemical reinforcement and the particle-contact reinforcement. The physical explanations of these mechanisms are presented as follows.

Volume-filling reinforcement
The volume-filling reinforcement can be explained as the stiffening due to the disturbance of the stress/strain fields in the soft matrix causing by the addition of stiff inclusion particles [86].
When the particles' concentration is very low, the disturbed area caused by each particle does not interact with each other, Fig. 13a; while with the concentration increasing, the disturbed areas caused by different particles may overlap and interact with each other, which is called as the ''inter-particle interaction" as mentioned in the previous sections, Fig. 13b.
According to the above definition of the volume-filling reinforcement, it is obvious that the stiffening effect of this mechanism is dependent on the volume fraction of the particles. In addition, the geometrical properties of the particles (i.e. the size, the shape, the angularity, etc.) make a major contribution to the stiffening effect of the volume-filling reinforcement as well.

Physiochemical reinforcement
The physiochemical reinforcement is defined as the stiffening because of the physicochemical interactions (i.e. absorption, adsorption, etc. [71]) at the interface between the matrix and inclusion particles. These interactions yield coating layers around the inclusion particles which increase the composite's stiffness [87], Fig. 14.
The stiffening effect of the physiochemical reinforcement is mainly affected by the geometrical and mineral characteristics of the inclusions. High surface areas, rough surface textures and high surface activities contribute to the increase of the composite's stiffness [88][89][90][91][92].

Particle-contact reinforcement
The particle-contact reinforcement refers to the stiffening resulting from the contacts among different particles [68,71]. When the concentration of particles is low, the particles are randomly distributed within the matrix and do not contact each other, Fig. 15a. Whereas, with the increase of the particles' concentration, a group of particles start touching each other and gradually form a skeleton framework [68,86], see Fig. 15b. Due to the formation of the skeleton framework, the stiffness of the composite becomes much higher than the bulk matrix.
It is obvious that the stiffening effect of the particle-contact reinforcement depends on the particles' concentration. Apart from that, it also depends on the loading condition, the temperature/frequency of the material, the geometrical properties of the particles, etc. For example, the particle-contact reinforcement is supposed to be more pronounced in a composite consisting of more angular particles with rough textures.
Overall, three different mechanisms result in the stiffening of an asphalt material. It is highlighted here that at a certain condition, it is possible that all of these mechanisms simultaneously play important roles in the stiffening of the material. It is also possible that the material's overall behavior is dominated by only one mechanism or two mechanisms, and the stiffening effects of the other(s) can be neglected. Therefore, in order to effectively predict an asphalt material's properties, it is necessary to figure out the dominant stiffening mechanism(s) beforehand.

Possible explanation of the models' limitations
On the basis of the understanding of the stiffening mechanisms of asphalt materials, the limitations of the commonly used CBMM at high concentrations and high temperatures can be explained as follows: None of these models can explicitly account for the interparticle interaction at high concentrations [86]. The Dilute model just describes the stiffening effects due to one particle, and thus, there is no interaction considered. Other models bring in the inter-particle interaction. However, since these models just take the set of all the particles as one phase, the particles' locations or their relative configurations are not taken into account in the predictions. In addition, some factors that impact the stress and strain distributions, such as the particles' size, irregular shape, angularity, etc., are not considered as well. There are no physicochemical interactions and particle contacts considered in these models; therefore, neither the stiffening effect of the physicochemical reinforcement nor that of the particle-contact reinforcement can be captured by these models.

Efforts on improving the accuracy of the upscaling results
As can be seen from the above section, it is necessary to take different stiffening effects into consideration so as to improve the predictions' accuracy. For this purpose, researchers have made efforts to either modify the commonly used CBMM or develop new models. In this section, some of these modified and new developed models are discussed.

Models with the consideration of aggregates' size
In asphalt mixture, different grades of aggregate particles are used. It has been generally found that the size of aggregate particles affects the mechanical properties of the mixture [93]. However, in commonly used CBMM, only the inclusions' volume fraction is considered; and there are no parameters representing the properties of inclusions in the length scale. In order to take the effect of the aggregates' size into account, Li et al. [94] have developed a three-phase micromechanical model (named as Li's model in this study).

Li's model
Similar to the GSC model, Li's model assumes a matrix-coated circular inclusion to be embedded into an equivalent composite medium, see Fig. 16. However, unlike the GSC model, the equivalent medium in Li's model is finite rather than infinite. The other main difference between these two models is the ratio of the inclusion's radius to that of the matrix. In the GSC model, this ratio is a. A low volume fraction of particles b. A high volume fraction of particles assumed to be only related to the volume fraction of the inclusion particles, while in Li's model, it is also dependent on the size and the gradation of these particles. The effective Young's modulus, E eff , can be derived by using the elasticity theory with radial stress uniformly applied at the boundary, see Equation (81). The value of E eff is a function of the mechanical properties of each phase, the inclusion's size and the matrix's thickness.
where E, m and R denote the Young's modulus, the Poisson's ratio and the radius, respectively; the subscripts ''m" and ''i" represent the matrix and the inclusion, respectively; and m eff is the Poisson's ratio of the mixture.
The values of E i , E m , m i and m m are known once the materials are selected. The value of m eff can be estimated according to each phase's volume fraction and mechanical properties [95]. To determine the values of R i and R m , it is assumed that the matrix coats each particle with the same thickness. With this assumption, the value of R m -R i can be determined from the phases' volume fractions and the gradation of aggregates, Equation (82).
where n is the total grades of aggregates by sieving, i.e. 0.075 mm, 0.15 mm, 0.3 mm, 0.6 mm, etc.; R j is the average radius of aggregates between the j th grade and the (j + 1) th grade; and P j and P j+1 are the weight percentages of aggregates passing through the j th grade and the (j + 1) th grade, respectively.
By combing Equation (81) and (82), a can be rewritten as Equation (83). It can be seen that the value of a is dependent on the gradation of aggregates. For different sizes of aggregates, the corresponding values of predicted E eff (R j ) are also different. Each value of E eff (R j ) can be explained as the contribution made by a particular size of aggregate particles. By adding the contributions from different sizes of aggregates, the overall effective modulus can be obtained, see Equation (84).
Li's model has been further improved by Shu, Huang et al. [50,96,97] to a four-phase model by adding another layer between the inclusion and the matrix. They have also improved Li's model from a two-dimension model to a more reasonable threedimension model. By using the same approach of considering aggregate gradation, they have taken the effects of the size distribution of air voids into consideration as well.

Performance of Li's model considering aggregates' size
Using Li's model, researchers [50] observed that the predicted modulus increased by around 20% when the maximum aggregate size increased from 4.75 mm to 19 mm. However, this limited increase of modulus was not enough to account for the significant difference between the predictions and the experimental results at lower frequencies (i.e. the measured modulus was around 100 times higher than the predicted value at a low frequency of 10 À3 Hz). This indicates that the effect of the aggregates' size may not take a leading role in the stiffening of asphalt mixture at low frequencies.

Models with the consideration of particles' configurations
To consider the inter-particle interaction more accurately, researchers [86,98,99] have proposed different models to bring in the effect of particles' configurations. In this section, two of these models which have been used for asphalt materials are introduced.

Physical interaction model
The physical interaction model [86] simulates the microstructure of a composite by idealizing the composite as a 2-D material, see Fig. 17. The particles are assumed as circles, and their size distribution is introduced into the model according to the gradation of the particles in the composite.
The stress disturbance function, which describes the stress state induced in an infinite matrix with a circular inclusion by remote stress, is given as Equation (85). It can be seen that the stress distribution around the particle is dependent on the inclusion's and the matrix's properties as well as the distance from the inclusion. An example of a typical r-r f curve is shown in Fig. 18. When the value of r approaches to 1 (a position close to the edge of the particle), the value of r f is quite high; while when the value of r is much higher than 1 (a position far away from the particle), the value of r f is close to the remote stress. Fig. 17. Illustration of the physical interaction model.

Finite effective medium
Inclusion Matrix R i R m Fig. 16. Illustration of Li's model.
where G 1 , G 2 , m 1 and m 2 are described in Equation (36); r is the normalized radial coordinate with respect to the radius of the inclusion; r f , which is called stress factor, is the ratio of the remote stress to the stress at the position of r.
The mixture rule which relates the composite's modulus to each phase's properties is illustrated in Fig. 19. Overall, the inclusion phase and the matrix phase are arranged in series, which indicates that the average stresses of these two phases are the same. Therefore, the value of G eff can be given as Equation (86).
G eff1 in Equation (86) is the effective shear modulus of the matrix by considering the stiffening effect of different inclusion particles on the individual volumes of the matrix. The moduli of these individual volumes are determined according to the stress distributions in the matrix. The stress factor for the matrix's volume surrounding a given aggregate particle r f is given by averaging the stress between that particle and each of its nearest neighbors, see Equation (87). The stress factor between two particles is calculated by taking the product of the individual stress factors for the two particles at their separation midpoint. The number N denotes the number of surrounding particles. Once the stress is computed for each particle, the modulus of the individual volumes is given as the weighted average of the stress factors, see Equation (88).
where (r f ) ij and (r f ) ji denote the stress factors at the midpoint between particle i and j introduced by particle i and particle j, respectively; and r ij and r ji are the normalized radial coordinates of the midpoint with respect to the radius of particle i and particle j, respectively.
where A k is a value of r f ; P k is the proportion of the total matrix volume with a stress factor of A k ; M is the number of different values of A k .

Ju and Chen model
Ju and Chen [98] have proposed a micromechanical model (termed as ''J-C model") to predict the effective moduli of a two-phase composite containing randomly distributed spherical particles by considering the interaction between two particles x 1 and x 2 , see Fig. 20a.
In order to determine the effect of the interaction, it is required to know the distance between two particles. However, for a composite with randomly dispersed particles, it is difficult to determine the exact distance between any two particles. Therefore, the conditional probability function P(x 2 |x 1 ), which describes the probability of finding x 2 at a given distance r, is introduced, see Fig. 20b. The value of P(x 2 |x 1 ) can be given as Equation (89).
where a denotes the particles' radius; r denotes the distance between two particles' centers; N is the total number of the particles in the composite with a volume of V; and g(r) is known as the radial distribution function, which describes the configuration of the particles. For different radial distribution functions, the values of g(r) are different. For example, for a statistically uniform radial distribution function, g(r) is given as one, while for the Percus-Yevick (P-Y) radial distribution function, the value of g(r) is a function of the particles' volume fraction.
With the consideration of the pairwise particle interaction and the conditional probability function, the values of K eff and G eff by using the J-C model are finally given as:

Performance of models considering particles' configurations
In comparison to the commonly used CBMM, both the physical interaction model and the J-C model enable the consideration of the particles' configuration and the inter-particle interaction. Furthermore, the physical interaction model includes the effect of different particles' sizes.
Researchers [86] applied the physical interaction model to two types of mastic with various filler concentrations. It was shown that the commonly used CBMM provided an accuracy of a filler concentration of around 10%. By contrast, the physical interaction model performed quite well in predicting the magnitude of the shear modulus up to a filler concentration of approximately 40%. Beyond this point, the model began to underestimate the measured modulus since it cannot account for the stiffening effect of the particles' contacts.
The J-C model was utilized to estimate the modulus of dense asphalt mixture by upscaling from the mortar scale [100]. The authors observed that the predicted moduli using the J-C model with the P-Y radial distribution function were much higher than those predicted using the MT model and the Differential model. For example, at lower frequencies, the J-C model provided a modulus of almost three times higher than the other two models. Compared to the experimental results, the J-C model predictions were found to be accurate at both high frequencies and low frequencies.
It is noted that the good performance of the J-C model may be due to the fact that there are not a large number of coarse particles' contacts formed in the dense asphalt mixture. However, this may not be the case for other types of asphalt mixture, such as porous asphalt mixture and stone matrix asphalt, where the volume fractions of course aggregates are much higher. Therefore, the applicability of this model on other types of mixtures requires to be further validated. In addition, the J-C model only provides the analytical solution for a two-phase composite. The use of this model to a multi-phase composite requires the sequence step method. The limitation of this method, as mentioned earlier, is that the predictions depend on the order of adding different phases.

Models with the consideration of physicochemical interactions
As mentioned earlier, due to the physicochemical interaction, there exist coating layers forming at the surface of the inclusion particles. In order to take the physicochemical interaction into account, researchers [70,71] have modified the commonly used CBMM by adding one layer between the inclusion and the matrix. Depending on the assumed different properties of this layer, the modified GSC model and a four-phase model have been proposed. In the modified GSC model [71], the modulus of the coating layer is considered to be the same as the inclusion's modulus. This indicates that the coating layer acts as an extension of the inclusion, see Fig. 21. Therefore, the effective volume fraction of the inclusion / e can be calculated by adding the real volume fraction of the inclusion / 2 and the volume fraction of the coating layer / c , see Equation (92).
It is assumed that the ratio of / c to / 2 (represented as a) is independent of the value of / 2 , see Equation (93). Therefore, the value of a can be calibrated from a low concentrated condition when the effect of the particle-contact reinforcement is not significant. Once a is known, / e at higher concentrations can be obtained using Equation (94).

5.3.2.
A four-phase model with a non-rigid coating layer 5.3.2.1. Description of the four-phase model. In a four-phase model [70] of considering the effect of the physicochemical reinforcement, see Fig. 22, the shear modulus of the coating layer G c is assumed to be: (1) affected by temperature and frequency but independent of the coating material's volume; and (2) a value between the modulus of the original matrix G 1 and the modulus of the inclusion particles. On the contrary, the modulus of the residual matrix G 1r is assumed to be (1) affected by the frequency and temperature as well as the volume of the coating material; and (2) a value less than G 1 . The relationship between G c , G 1r and G 1 is given as: where / c and / 1r are the volume fractions of the coating layer and residual matrix, respectively. The thickness of the coating layer is assumed to be related to the volume of the matrix. The so-called ''maximum potential film thickness" [70] can be achieved when the matrix's volume is sufficiently high, see Fig. 23a. When the volume fraction of the matrix is low, there are not enough components to be adsorbed by the inclusion particles, which results in the reduction of the thickness of the coating layer, see Fig. 23b.

Performance of models considering physicochemical interactions
The above descriptions show that both the modified GSC model with a rigid coating layer and the four-phase model with a nonrigid coating layer can account for the stiffening effect due to the physicochemical interaction. The modified GSC model is easy in implementation. However, its assumption of inclusion-like mechanical properties and a constant value of a may not be realistic. Also, the change in the properties of the residual matrix due to the formation of the coating layer is not considered in the model.
On the other hand, the assumed non-rigid coating layer and the changeable residual matrix in the four-phase model are more reasonable and realistic. However, in this model, there are many unknown parameters to be determined, i.e. the maximum potential film thickness, the modulus of the coating layer, etc. Since it is very difficult to obtain these parameters experimentally or theoretically, they have to be adjusted by conducting a series of complex procedures during the process of the upscaling.
Researchers [70,71] observed that comparing against the commonly used CBMM which provided an accuracy of around only 10%, the accuracy of the predictions improved with the consideration of physicochemical interactions. In their research work, authors also highlighted that there was an optimal value of the filler concentration (approximately 40%), beyond which the accuracy of predictions significantly decreased. This may indicate that these two models are still not powerful enough to represent the effect of the physicochemical reinforcement. This may also indicate that the effect of the physicochemical reinforcement, to a certain extent, accounts for the stiffening of the mastic, but it does not play the main role when the volume fraction and the test temperature are quite high.

Models with the consideration of particles' contacts
In order to take the particles' contacts into account, one of the commonly used methods is to modify the commonly used CBMM by combing them with the percolation theory [68,[101][102][103][104][105]. Another one is to quantify the internal structuralization that exists within a composite via laboratory tests and proceed to build the relationship between the internal structuralization and the modulus of the composite [106]. These two methods are discussed as follows.
5.4.1. Models using the percolation theory 5.4.1.1. Description of models using the percolation theory. The percolation theory [107,108] accounts for the effect of the particlecontact reinforcement by considering the non-percolated matrix within the clusters of connected particles, see Fig. 24.
By using the percolation theory, researchers [101,102] have constructed a four-phase model where the non-percolated matrix, surrounded by the inclusion phase and then the percolated matrix, is embedded into an infinite effective medium, see Fig. 25. The volume fraction of the non-percolated matrix / npm can be calculated by using Equation (96).
where / i and / m represent the volume fractions of the inclusions and the matrix, respectively; / a , known as the percolation threshold, is a critical value of / i when the non-percolated matrix begins to occur; and / i;max is the maximum value of / i beyond which no more inclusion particles can be packed into a given volume of the matrix.
Other researchers [68] have simplified the four-phase model to the modified GSC model (Fig. 26) by assuming the non-percolated matrix as part of the inclusion. In this case, the values of / npm and / i can be summed up as the effective volume fraction of the inclusion / eff . Furthermore, in some studies, the similar concept of considering particles' contacts as the percolation theory is used. For example, Lewis and Nelsen [109][110][111] generalized the Kerner equation [47] by introducing / i;max . Lackner et al. [49,80,112]  In the microstructural association model, the effect of the particle-contact reinforcement is accounted for by the internal structuralization of a composite [106]. The degree of the internal structuralization is qualified by using the structuralization index (SI). The value of SI can be calculated on the basis of the packing properties of the raw aggregates in the composite, see Equations (97) and (98).
where %Voids loose and %Voids comp are the contents of the voids in the loose and the compacted aggregates, respectively, which can be obtained from independent laboratory tests [113][114][115]; % Voids structure represents a critical value of the voids content in the packing aggregates when a stable structure begins to form; %Voids refers to the voids content in the packing aggregates of the given composite; and the parameter B describes the compaction needed to reach a stable structure, which requires to be calibrated from laboratory tests.
where k is the proportion of the aggregates that are in contact with each other. In Equation (99), G SM represents the contribution of the particle-contact reinforcement. It is assumed that the value of G SM depends on temperature and frequency which can be calculated by using a logistic function, Equation (100).
where a, b, v and d are fitted parameters, and G 1 is the matrix's shear modulus.

5.4.2.3.
Relationship between SI and the modulus of asphalt mixtures. The establishment of the relationship between SI and a mixture's modulus is on the basis of the finding that a consistent bilinear relationship exists between the value of SI and the stiffening ratio (SR) of asphalt materials at different length scales [70]. This indicates that the SR-SI relationship developed at a lower scale can be used to obtain the stiffening ratio of materials at higher scales. Therefore, the SR-SI relationship at the mastic scale can be first developed. Based on the experimental results of G eff of mastic, the values of G SM are calculated using Equation (99). By fitting G SM using Equation (100), the values of a, b, v and d can be determined. The values of a, b, v and d are fitted as functions of SI, which are further used to predict the modulus of asphalt mixture.

Performance of models considering particles' contacts
Researchers [68] evaluated the percolation theory by comparing the predicted mastic modulus using the modified GSC model with the non-percolated matrix layer to that using the original GSC model. The accuracy of the predictions using the original GSC model was found to be limited to a filler volume concentration of approximately 10%. By contrast, using the percolation theory, the accuracy could improve up to a concentration of 30% in which the values of / i;max and / a were adjusted. However, some of the calibrated values of / i;max were not in agreement with those measured from experiments; and the calibrated values of / a did not seem to be realistic. Other researchers [70] used the percolation theory in the four-phase model (see the section of ''A four-phase model with a non-rigid coating layer") attempting to account for the stiffening mechanisms occurring at high filler concentrations (>40%). Although a better match between the predictions and the experimental values was observed at frequencies higher than 10 À2 Hz, the model still failed to accurately predict the modulus at lower frequencies.

a. High volume of matrix b. Low volume of matrix
These undesirable predicted results indicate that the use of the percolation theory may not be able to account for the effect of the particle-contact reinforcement. This may be due to the fact that these models simply treat the connected particles as an entire phase while they never directly consider the contacts among different particles.
Comparing to the percolation theory, the use of SI in the microstructural association model is a more reasonable approach for capturing the particle-contact reinforcement in a composite. Researchers [106] showed that the mastic modulus at high filler concentrations (>40%) could be accurately fitted in the whole frequency range by combining Equation (99) and Equation (100). Furthermore, using the fitted relationships between the values of a, b, v and d and the value of SI in the mastic scale, the modulus of asphalt mixtures in a wide temperature range (from 20°C to 54°C) could be well predicted.
It is a novel work that the microstructural association model can provide accurate predictions at high inclusion concentrations and high test temperatures. However, the establishment of all the equations in the model requires a large number of laboratory tests, i.e. the shear modulus of mastic at different concentrations, the voids contents in the loose and compacted aggregates, etc. In addition, in the model there are many parameters required to be back-calculated, i.e. the value of B in Equation (97), the value of k at each concentration in Equation (99), the values of all the parameters in the four-phase model and the values of a, b, v and d in Equation (100). Overall, the implementation of the microstructural association model is tedious and complicated.

A summary of the modified and new-developed micromechanical models
A summary of the above discussed modified and newdeveloped micromechanical models is shown in Table 2. It can be concluded that the models accounting for the volume-filling reinforcement and the physicochemical reinforcement fail to improve the predictions at high concentrations and high temperatures significantly. This indicates that the behavior of highly concentrated asphalt materials at high temperatures/low frequencies may be mainly dominated by the particle-contact reinforcement.
The dominance of the particle-contact reinforcement at high concentrations and high temperatures could be explained on the basis of the hypothesis that in a highly concentrated asphalt composite, particles are close to each other, and the matrix layer located among them tends to be quite thin. When an external load is applied at lower temperatures, the matrix layers are stiff enough to avoid two adjacent particles moving closer to form direct contacts, see Fig. 27a. However, at high temperatures, since the matrix layers are soft, it would be much easier for the particles to move closer and form a considerable number of direct contacts, see Fig. 27b. Compared to the soft matrix layers, these direct contacts would play a leading role in transferring the load among different particles. Therefore, at high temperatures, the overall behavior of a highly concentrated composite is governed by the particle-contact reinforcement. Without considering this reinforcement, it is impossible for any CBMM to further improve the accuracy of the predictions.
The failure of the modified CBMM using percolation theory indicates that these models, by nature, cannot capture the stiffening effect of the particle-contact reinforcement [71,116]. This is due to the fact that in such models, the set of all the individual particles is simply represented as one inclusion. In this case, it is not possible to consider any characteristics of individual particles and the contacts between them. Furthermore, the use of SI in the microstructural association model provides a way to capture the properties of contacting aggregates and their effects on the behavior of asphalt mixtures. However, since SI is obtained from laboratory tests, it cannot give any fundamental insight into the physical mechanisms related to the particle-contact reinforcement.

Research challenges and future recommendations
From the above discussion, it can be concluded that the most challenging work about the application of CBMM for asphalt materials is to capture the effect of the particle-contact reinforcement at high concentrations and high temperatures. For this purpose, the authors proposed the following approaches.

Reasonable and realistic microstructure modeling
Particle-contact reinforcement is related to the microstructure of asphalt materials. A reasonable and realistic microstructure is a precondition for accurately capturing the characteristics of particles' contacts.
For example, the microstructure modeled by the physical interaction model can be considered as a reasonable and realistic microstructure for asphalt mastic and dense asphalt mixture where the inclusions are surrounded by a continuous matrix. However, it cannot be considered as a suitable microstructure for porous asphalt mixture where the matrix tends to discontinuously locate between different particles [17]. In this case, a discretebased micromechanical model where a coated particle-onparticle skeleton is bonded by discontinuous matrix layers may be more reasonable and realistic, see Fig. 28.
Various researchers have proposed different approaches to estimate the mechanical properties of a composite using discretebased micromechanical models [117][118][119][120]. Although the applications of these models on bonded granular materials such as glass beads packs, frozen sand, etc. have been investigated [117,[120][121][122], limited research work has been conducted to use these models for asphalt materials. Relevant studies can be found in the work of Cheung et al. [123], Zhu and Nodes [122], etc., where a discretebased micromechanical model was used to simulate the creep characteristics of asphalt mixtures. However, these studies did not provide a rigorous way to implement the model for asphalt materials especially in terms of the determination of each phase's geometric characteristics, such as the radius of the particles, the radius and thickness of the binder layer, etc. Therefore, more efforts are required to investigate the implementation and evaluation of these models for asphalt materials.

Consideration of physical mechanisms behind the particle-contact stiffening effect
In order to describe the particle-contact reinforcement, it is necessary to understand the physical mechanisms related to the contacts of the particles, such as friction, interlocking, Hertzian normal contact, etc. The properties of these mechanisms are associated with not only the properties of the composite such as the contents of each phase, the gradation of aggregates, the morphology of aggregates, etc. but also the loading conditions such as compressive or tensile loads, with or without confinement, test temperatures/frequencies, etc.
For example, for asphalt mixtures consisting of angular aggregate particles with a large number of textures, the friction and interlocking mechanism may make considerable contributions to the behavior of the mixture. However, when the mixture contains round and smooth particles, the contribution from the friction and the interlocking may be neglected. Physical mechanisms related to particle contacts have been widely realized by researchers [124,125]. However, the challenge is to quantitatively evaluate their effects on the mechanical properties of asphalt materials using an analytical way. Currently, most studies in the pavement field addressed this issue by just using calibration factors [125,126]. Since calibration methods do not give any fundamental insight into how these mechanisms affect the characteristics of particles' contacts, a more logical and reasonable method needs to be developed in future research studies.  Although characteristics of particles' contacts have not been thoroughly studied in the pavement filed, a large amount of relevant research work has been conducted in other fields, such as concrete, granular solids, etc. Researchers [127,128] have developed different analytical models to quantify the forces transferred through Hertzian contacts for granular solids. These developments could be adapted for asphalt materials in the near future.

Consideration of matrix's viscoelastic behavior
As mentioned earlier, the micromechanical models which are originally developed for elastic materials can be used for viscoelastic materials according to the elastic-viscoelastic correspondence principle. However, the viscoelastic properties of the matrix may lead to the change of the connected particles' behavior. In this case, the way of simply replacing the elastic modulus of the matrix with a complex number can never capture this change.
For example, when a cyclic compressive force is applied to an asphalt mixture at a high temperature, initially, the load is mainly transferred via the binder layer. However, due to the viscoelasticity nature of the binder layer, the displacement of the binder layer increases with the increase of the loading time. With the increase of the displacement of the binder layer, the distance between adjacent particles decreases. At a critical time, many particles start to contact each other, and thus the particles' direct contacts start to transfer the load. In order to describe the change of the load transfer mechanism among different particles, the viscoelastic behavior of the binder layer needs to be considered in the model to describe its time-and force-dependent displacement.
To date, there are few analytical micromechanical models considering the effect of the viscoelastic properties of the binder layer on the characteristics of particle contacts in asphalt materials. Nevertheless, there have been attempts to model the effect of the viscoelastic properties of the binder layers on the evolution of the microstructure of an asphalt mixture [123,129]. In the proposed models, the mixture was represented as a set of discrete spherical particles bonded by thin binder layers, as shown in Fig. 28. With the assumption that the particle centers moved compatibly with the macroscopic strain rate, the microscopic deformation of a binder layer was estimated. Combining these studies with a model describing the behavior of contacting particles (as discussed in the previous subsection), it is supposed that the evaluation of the effect of the viscoelastic properties of the binder layers on the characteristics of particles' contacts can be achieved.

Utilization of DEM
On the basis of the above discussion, it can be concluded that to accurately capture the stiffening effect of the particle-contact reinforcement, it is required to model the microstructure of a material in a reasonable and realistic way. In addition, it is required to accurately describe the interactions between any two particles under a given loading condition. In order to achieve this, the morphological characteristics (i.e. shape, size, orientation, texture and angularity) and the spatial distribution of each particle need to be accurately represented. However, it is almost impossible to describe these complex characteristics by using an analytical micromechanical model. Most of the existing micromechanical models can only consider the shape of particles as regular shapes, i.e. spheres, ellipsoid or cylinder, and it is almost impossible to bring in the angularity and the texture of the particles into these models. Therefore, an accurate prediction of the particle-contact reinforcement requires a robust numerical model. Based on the above realization, it is considered that DEM is a promising technique for accurately capturing the stiffening effect of the particle-contact reinforcement. This is because in DEM, var-ious material phases can be modeled as assemblies of very small discrete elements. By modeling individual particles as a clump of small discrete elements, it is possible to model complex morphology of aggregates [7,130]. Furthermore, different contact constitutive models can be used in DEM to describe the interactions between any two discrete elements [7]. By combing the characteristics of each discrete element and their interactions, the overall behavior of asphalt materials under different complex conditions is supposed to be captured.
Over recent decades, the application of the DEM technique to simulate the mechanical behavior of asphalt materials has been widely studied [7,[131][132][133][134][135][136][137]. The generation of a realistic microstructure and the utilization of different contact constitutive models can be found in many research studies [130,132,133,[136][137][138]. Since the review of the DEM technique is beyond the scope of this paper, details about these studies will not be further described.
Overall, four recommendations from different perspectives have been made to accurately capture the stiffening effect of the particle-contact reinforcement. A summary of these recommendations is given in Table 3.

Conclusions
This study provided a comprehensive review of CBMM (the first part) and their applications for asphalt materials (the second part). In the first part, the commonly used CBMM were introduced in detail and their advantages and disadvantages were discussed and compared. Since CBMM are initially developed for elastic composites, a general procedure to apply them for viscoelastic composites was elaborated. The main findings of this part include: Mainly two different approaches are used to estimate the effective properties of a composite: (1) Geometry based approach to obtain a closed-form solution for the effective properties; (2) Bound based approach to determine the lower and upper bounds of the possible range of these properties. Models (i.e. the Voigt model, the Reuss model, the Hirsch model, etc.) developed from an arrangement of individual phases in parallel or series or a combination of them can provide a closed-form solution for the effective properties. However, the calibration of the proportion of the series and the parallel part of each phase is always required for accurate predictions. Table 3 A summary of future recommendations.

Recommendations Reasons
Reasonable and realistic microstructure modeling A reasonable and realistic microstructure is a precondition for accurately capturing the characteristics of particles' contacts. Consider physical mechanisms behind the particle-contact stiffening effect Depending on composites' characteristics and loading conditions, different mechanisms are related to particles' contacts. Consider matrix's viscoelasticity Viscoelastic properties of the matrix may lead to the change of the connected particles' behavior.

Utilize DEM
The complex morphological characteristics and the spatial distribution of aggregate particles can be described using DEM; Different contact constitutive models in DEM can describe the interactions between any two discrete elements.
Eshelby-based micromechanical models are developed based on the Eshelby's solution where an ellipsoid inclusion is embedded into an infinite matrix. The Dilute model neglects the interparticle interaction while other models (the SC model, the GSC model, the (N + 1)-phase model, the MT model, the Differential model and the Composite Sphere model) consider the interaction by using different ways. The H-S bounds can provide better estimates for the effective moduli of a composite in comparison to the Paul's bounds. However, neither of them provides narrow bounds results when the differences of moduli between individual phases are significant. CBMM can be used for viscoelastic composites according to the elastic-viscoelastic correspondence principle. In the frequency domain, the complex moduli of a viscoelastic phase can directly replace the corresponding elastic moduli. In time domain, the relaxation moduli or creep compliances of a viscoelastic phase require being transformed into the LC domain to calculate the effective relaxation moduli or effective creep compliances.
In the second part, a concise summary was given about the application of commonly used CBMM for asphalt materials. In order to explain the limitations of these models, the stiffening mechanisms of asphalt materials were introduced. Furthermore, the modified and new developed CBMM for improving the accuracy of the predictions were summarized and discussed. In the end, research challenges were highlighted and further recommendations were proposed. Overall, the obtained conclusions are given as follows: The performance of the models is dependent on the volume concentration of the inclusions and the test temperature. At very low concentrations, most of the models can obtain good predictions. At high concentrations and low temperatures, the predicted results differ between different models. At high concentrations and high temperatures, none of these models can obtain good predictions, and in general, the predicted modulus is much lower than the measured value. The upscaling from the properties of a higher scale matrix is more accurate than upscaling from a lower scale matrix. The multiscale modeling technique can increase the predictions' accuracy as well, but the improvement is not very significant at high concentrations and high temperatures. Generally, three mechanisms contribute to the stiffening of asphalt materials: the volume-filling reinforcement, the physiochemical reinforcement and the particle-contact reinforcement. The commonly used CBMM consider but not explicitly account for the effect of the volume-filling reinforcement, while the effects of other mechanisms are not even considered by these models. The consideration of the effects of the aggregates' size, the particles' configuration and the coating interface cannot significantly improve the accuracy of the predictions at high concentrations and high temperatures. This indicates that the behavior of highly concentrated asphalt materials at high temperatures may be mainly dominated by the particle-contact reinforcement. The modified CBMM using the percolation theory cannot account for the particles' contacts because these models simply treat the connected particles as an entire phase while they never directly consider the contacts among different particles. On the other hand, the use of SI in the microstructural association model is a more reasonable approach for capturing the particle-contact reinforcement. However, the establishment of all the equations in the model requires a large number of laboratory tests; meanwhile, most of the parameters in the model have to be calibrated. The most challenging work about the application of CBMM for asphalt materials is to capture the effect of the particlecontact reinforcement at high concentrations and high temperatures. For this purpose, the authors recommended building a realistic and reasonable microstructure in the models. Moreover, it is crucial to figure out the mechanisms behind the stiffening effect of particles' contacts and to account for the matrix's viscoelastic properties. However, considering the limitations of analytical models in describing the complex morphological characteristics of individual particles and the interactions between any two particles, it was further recommended to use discrete element model to predict the behavior of asphalt mixtures under different complex conditions.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.