Contact mechanics based solution to predict modulus of asphalt materials with high porosities

Article history: Received 19 December 2020 Revised 11 April 2021 Accepted 17 April 2021 Available online 22 April 2021


Introduction
Due to the advantage of improving wet skid resistance and reducing hydroplaning and ''splash and spray" effects, asphalt mixtures with high porosities, known as porous asphalt (PA) mixes, have attracted the attention of many researchers and engineers [1][2][3][4]. PA mixes can be considered as environmental-friendly materials since they are capable of reducing pavement-tire interaction induced noises. All these advantages motivate the wide use of PA mixes in some countries such as the Netherlands where PA mixes surface more than 80% of the highways [5].
In order to come up with a proper PA mix design, a good estimation of the mechanical properties of a given mix is often required. Commonly available pavement analysis tools such as 3D-move [6], together with the mechanical properties of the mix, can be used to estimate the lifetime expectancy of the designed pavements under critical loading conditions. In practice, in order to determine a mix's mechanical properties, a series of laboratory tests are often conducted. However, such laboratory tests are usually not only expensive but also time-consuming.
In order to reduce experimental work, researchers [7,8] have made efforts to obtain mechanical properties of an asphalt mixture directly from the properties of its constituents. These empirical/ semi-empirical relationships have been broadly used for estimating the mechanical properties of asphalt mixtures [9,10]. However, such models contain various parameters that are calibrated on the basis of a large amount of experimental data. The need for calibration limits the applicability of empirical/semi-empirical models as the determined values of the parameters are only applicable to the mixtures they are calibrated for [11].
With the above realization, researchers started to look for theoretical relationships to predict the mechanical properties of asphalt mixtures. In recent decades, micromechanical models which are developed on the basis of basic mechanical theories (such as continuum mechanics) have drawn the attention of many researchers [12][13][14][15][16][17][18]. Among different types of micromechanical models, numerical micromechanical models, which are developed based on finite element method (FEM) and/or discrete element method (DEM), have been favored by many researchers [12][13][14][15] mainly because these models can describe almost realistic microstructures of asphalt mixtures (i.e. different sizes and shapes of aggregates). However, the description of such complex microstructures requires the development of a huge number of finite or discrete elements. This indicates that powerful computational facilities are always needed to analyze the mechanical properties of asphalt mixtures, and it always takes a long time to obtain the predicted results.
Furthermore, although the microstructure of asphalt mixtures can be realistically described using FEM or DEM, accurate predictions still rely on the calibration of some input parameters. For example, experimental results are typically needed to calibrate the parameters of the contact constitutive models between two particles when DEM is used to simulate the mechanical properties of asphalt mixtures [12,13].
On the other hand, analytical micromechanical models can be a promising alternative to determine the mechanical properties of asphalt mixtures in a more effective way. These models simplify the descriptions of the microstructure of asphalt mixtures. For example, aggregates are typically modelled as identical spherical or ellipsoidal particles. However, these simplified descriptions allow obtaining analytical solutions to the mixture's mechanical properties. Therefore, in comparison to FEM and DEM, analytical micromechanical models are much less time and facilitiesconsuming. Furthermore, although analytical micromechanical models cannot realistically describe the microstructure of asphalt mixtures, their capabilities of accurately predicting the mechanical properties of the mixtures have been reported in many research studies [19,20].
Currently, research work regarding analytical micromechanical models has mainly focused on the continuum-based micromechanical models (CBMM), such as the Self-consistent (SC) model [21], the Mori-Tanaka (MT) model [22], the Differential model [23], the generalized Self-consistent (GSC) model [24], etc. These models are developed based upon a matrix-inclusion system where an inclusion is embedded into a continuous matrix [25]. Although CBMM could provide accurate predictions for materials with low concentrations of particles, the predictions at high concentrations did not match well with the experimental results, especially at low frequencies [16][17][18].
Researchers postulated that the inaccuracies of the predictions could be related to the fact that CBMM do not sufficiently account for the stiffening effect of particles' interactions [16,17]. In a composite, particles' interactions are always associated with the characteristics of individual particles, such as their sizes, shapes, locations, etc. On the contrary in CBMM, all the particles are represented as one phase, which makes it impossible for the models to include these characteristics.
In recent years, few micromechanical models where a specific microstructure of a composite can be defined have also been used by pavement researchers [26]. For example, in the physical interaction model [26], different sizes of individual particles are embedded into a continuous matrix, see Fig. 1. On the basis of the distribution of the particles and their different sizes, the stiffness of the composite can be estimated by quantifying the stiffening effects of individual particles and their interactions.
Although the physical interaction model provides a way to consider the characteristics of individual particles, it is only suitable for composites without or with few air voids, such as mastic and dense asphalt mixtures. This is because the matrix in the model is considered to be continuous. However, in a PA mix, since there is a high content of air voids, the matrix phase is expected to be discontinuously located among densely packed aggregate particles. Therefore, in order to accurately predict the mix's mechanical properties, a different micromechanical model which can reasonably describe the microstructure of the mix is required.
In comparison to the physical interaction model, the discretebased micromechanical model (DBMM) [27][28][29] potentially provides a better way to describe the particles' interactions in PA mixes. This model simulates a bonded granular material as a particle-on-particle skeleton bonded by binder layers. This arrangement is naturally expected to represent a more reasonable and realistic microstructure for PA type mixes. Moreover, such a model has the capability to include several geometrical parameters to describe the characteristics of individual particles, for example, the radius of the particles, the average number of contacts per particle, etc.
Solutions to the effective moduli of a bonded granular material using DBMM have been derived by different researchers [27,29]. In this study, solutions provided by Dvorkin, et al [27] were used and thus DBMM will be referred to as Dvorkin's model in the following sections. Dvorkin's model has been successfully used to predict the modulus of granular materials bonded by elastic binders such as frozen sand, glass beads packs, etc. [30]. According to the elasticviscoelastic corresponding principle [31], models for elastic materials can be used for viscoelastic materials as well. Therefore, it is highly possible that Dvorkin's model could be used for estimating the modulus of asphalt materials. However, limited research studies have provided thorough investigations in this area. For example, Zhu and Nodes [28] used Dvorkin's model to simulate the creep behaviour of asphalt mixtures, limiting their study to some randomly assumed geometrical parameters, whereas, it is obvious that the values of geometrical parameters significantly affect the predictions.
In light of the above discussions, the main aim of this study is to investigate the performance of Dvorkin's model on predicting the modulus of PA mixes. The scope of this study includes: to propose a framework for predicting the modulus of PA mixes using Dvorkin's model; to evaluate the performance of Dvorkin's model based upon the experimental results of the complex modulus of PA mixes.
In order to make this paper understandable for readers, it has been arranged in different sections. The first and second sections will introduce and provide background information about Dvorkin's model. The proposed framework for predicting the effective moduli of PA mixes using Dvorkin's model will be elaborated in the subsequent section. In the fourth section, materials characteristics and performed laboratory tests will be presented. The fifth section will discuss the obtained model parameters. The last two sections will describe the predicted results and conclusions.

Dvorkin's model
In Dvorkin's model, the overall effective moduli of a granular material are derived by homogenizing the stiffness of an assembly of two-bonded particle systems that orient in different directions. This section first describes the stiffness of a two-bonded particle system, on the basis of which, a homogenization technique is introduced to obtain the effective (homogenized) moduli of a bonded granular material.
2.1.1. Stiffness of a two-bonded particle system A two-bonded particle system consisting of two equal-sized spheres and a layer of binder is illustrated in Fig. 2a. The spherical particles are assumed to have an arbitrary radius R, with the binder layer radius of a and a minimum height of 2 h 0 . It is noted here that generally, particles in a granular material have different sizes and irregular shapes. Since closed-form solutions for calculating the stiffness of a two-bonded particle system with irregular shapes and different sizes have not yet been accurately derived, past researchers have considered these particles as spheres with uniform sizes [27,29].
In Dvorkin's model, the stiffnesses of a two-bonded particle system S are defined as: where F is the applied force on the system; d is the displacement of the centre of one particle relative to the median plane (the x axis); and the subscripts n and s denote the normal and the tangential direction, respectively. The value of d is related to the displacement of the binder layer's surface (represented as V and U in the normal direction and the tangential direction, respectively) as well as that of the sphere's surface (represented as v and u in the normal direction and the tangential direction, respectively), see Eqs. (2) and (3).
In order to calculate the displacement of the binder layer, researchers [27] proposed a solution for a similar problem assuming a thin binder layer as an elastic base. In this solution, V(x) and U (x) were related to the stress along the interface r via Eqs. (4) and (5), respectively. Furthermore, since the contact region is small as compared to the particle, v(x) and u(x) can be approximated as the surface's displacement of an elastic half-space [27], see Eqs. (6) and (7), respectively.
VðxÞ hðxÞ ð4Þ are the Poisson's ratio and the shear modulus of the binder layer, respectively; m p and G p are the Poisson's ratio and the shear modulus of the particles, respectively; h(x) is the halfthickness of the binder layer, which can be approximated as By combining Eqs. (2), (4) and (6) or Eqs. (3), (5) and (7), the displacement of the binder layer's surface can be obtained giving a non-zero constant value of d, see Eqs. (9) and (10). These two equations are known as the Volterra integral equations of the second kind. Since this type of integral equations is too complicated to be solved analytically, a numerical technique was used, see Appendix A. Once the displacement of the binder layer's surface is known, the value of r can be computed from Eqs. (4) and (5). The value of F is determined by integrating r on the whole area of the interface, see Eqs. (11) and (12). Fig. 2. A two-bonded particle system.
It is noted that the above solutions were developed for particles bonded by an elastic binder layer. According to the elastic-viscoelastic corresponding principle, these solutions can be used for viscoelastic materials by replacing the elastic moduli of the binder layer with the complex moduli of the viscoelastic binder layer.

Homogenization technique
In order to upscale the mechanical properties of a two-bonded particle system to the effective moduli of a granular material, a homogenization technique is needed. As the first step, an assembly of bonded particles, as illustrated in Fig. 3, is considered as a representative volume element (RVE) of a bonded granular material. In the RVE, the particles are considered as spheres with a uniform radius of R; and all the binder layers are assumed to have a radius of a and a minimum height of 2 h 0 .
Two spheres r and s, which are initially in contact with each other, have the position vectors of X (r) and X (s) in the global coordinate system, respectively. Under a deformation u subjected to the boundary of the system, the centers of r and s also undergo displacements, represented as u (r) and u (s) , respectively. By using the kinematic hypothesis that the strains throughout the packing are uniform, the values of u (r) and u (s) can be calculated via Eq. (13).
where he ij i is the average strain applied on the packing. By symmetry, it can be obtained that the center of the median plane between the s th and the r th particles undergoes a displacement of (u (r) + u (s) )/2. Thus, relative to the median plane's center point, the displacement of the r th sphere center is (u (r) -u (s) )/2. This relative displacement can be separated into a normal component d n (Eq. (14)) and a shear component d s (Eq. (15)). The total force F on the s th sphere due to its contact with the r th sphere is given as the sum of the force in the normal direction and that in the tangential direction, see Eq. (16) Furthermore, it is defined that I is the unit vector along the line of centers of two particles: Substituting Eq. (18) into Eq. (17), the value of F can be further given as The average stress field hri of the granular material with a total volume V is related to the stresses within individual particles by using the following equation: where r ij (s) is the stress within the s th particle; V s represents the volume of the sphere and N p is the number of all the particles within V.
By using the divergence theorem, the value of r ij where S s is the surface of the s th particle; x' denotes the position vector of a point on S s relative to the center of the particle; s (s) denotes the traction across S s . From Fig. 3, it can be seen that the value of s (s) is nonzero only at the position where the particle contacts with other particles. Therefore, Eq. (21) can be further written as Z Vs r ðsÞ where n is the number of contacts per particle. Since the contact area between the s th and the r th particle is small, the value of x' at the contact area can be approximated by the position vector of the center of the contact area relative to the center of the s th particle: By substituting Eq. (23) into Eq. (20), the value of hri can be written as This value can be further expressed as Eq. (25)   Assuming that the geometry of the packing is statistically isotropic and the distribution probability of the contact points over the surface of a sphere is uniform, the summation in Eq. (25) can be represented in terms of averages: where the brackets <Á> denote the average over all uniformly distributed unit vector I. In Eq. (26), the total volume of the system V can be given as where V p and / p denote the volume and the volume fraction of the particles in the RVE, respectively. When all the particles are spheres with uniform sizes, the value of V p is given as Combining Eqs. (27) and (28), the value of N p /V can be expressed as Therefore, Eq. (26) can be further written as

Effective moduli of a bonded granular material
The stiffness tensor C* of a bonded granular material can be derived from Eq. (30) as According to the relationships in Eq. (32) and the definition of the stiffness tensor for an isotropic material in Eq. (33), the Young's modulus E* of a bonded granular material can be calculated using Eq. (34). By substituting the determined values of S n and S s from Eq. (1) into Eq. (34), the value of E* can be obtained.
where n denotes the coordination number; / p is the volume fraction of the particles; R represents the radius of the particles; h 0 refers to the minimum distance between adjacent particles; and S n and S s denote the stiffnesses of a two-bonded particle system in the normal and in the tangential direction, respectively.

Proposed framework to predict stiffness of PA mixes
The proposed framework for predicting the stiffness of PA mixes contains three main steps, see Fig. 4. A PA microstructure was assumed to be consisting of randomly packed spherical particles, mortar, and air voids. The total volume of mortar was categorized into two functional groups. Some mortar layers that are located between adjacent particles play a major role in binding particles together, and thus they were defined as the ''binding mortar layers"; whereas the remaining parts only play a role in coating individual particles, and thus they were defined as the ''coating mortar layers".
To predict the stiffness of PA mixes, the mechanical, volumetric and geometric properties of every phase are required. Generally, the mechanical and the volumetric properties can be directly measured from laboratory tests. Additionally, five geometric parameters need to be determined: (1) the radius of the spherical particles, (2) the thickness of the coating mortar layers, (3) the minimum thickness of the binding mortar layers, (4) the average coordination number, and (5) the radius of the binding mortar layers. These geometric parameters may also be measured using advanced technologies such as digital image processing. However, this method is difficult to be implemented since the values of all these geometric parameters vary with different locations. Thus, in this study, a different method that can be used to determine the geometric parameters of the proposed microstructure in a much easier but reasonable way was proposed, see Appendix B.
In the proposed microstructure, the properties of both the binding mortar layers and the coating mortar layers affect the stiffness of a PA mix. However, Dvorkin's model can only predict the stiffness of a skeleton framework consisting of particles, air voids and the binding mortar layers. Therefore, in the last step, a different methodology was employed to include the effect of the coating mortar layers in the model.

Proposed microstructure model for PA mixes
The microstructure model for PA mixes, see Fig. 5, was proposed according to the literature [32]. Individual particles (modelled as identical spheres with a uniform radius of R) are initially covered by mortar layers with a uniform thickness of t. When the mortar-covered aggregate particles pack together, the minimum distance between adjacent particles, represented as 2 h 0 , has an initial value of 2 t, see Fig. 5a. It was assumed that during the compaction process, the value of h 0 decreases, and depending on the amount of the compaction work, the value of h 0 ranges from t (without any compaction effort) to 0 (two particles contact each other), see Fig. 5b. It is noted here that during the compaction, not only the value of h 0 decreases but also the relative positions of aggregate particles change. However, since the change of the particles' positions is difficult to be analytically described, the proposed microstructure model did not take this phenomenon into account.
The volume of the binding mortar layer between two adjacent particles is illustrated in Fig. 6a. The binding mortar layers, with a radius of a, are assumed to be composed of two parts, see Fig. 6b. When two particles are compacted closer, their coating mortar layers overlap each other to form one part of the binding mortar layer with a radius of a 1 . Due to the overlap of the coating mortar layers, some mortar is squeezed out and forms the other part of the binding mortar layer with a radius of a 2 .

Procedures to calculate the stiffness of PA mixes
According to the literature [30], contributions to the stiffness of a mix made by the binding mortar layers and the coating mortar layers are different. The binding mortar layers play a major role in providing the stiffness of the mix because they enable the discrete particles to work together as a skeleton framework. By contrast, the coating mortar layers contribute to the stiffness of the mix by filling the air voids in the skeleton framework.
In light of the above discussions, the following procedure was proposed to calculate the stiffness of PA mixes. At first, the stiffness of the skeleton framework consisting of aggregate particles and the Step 1: Propose a microstructure for PA mixes (Section 3.1) Step 2: Determine the parameters in the proposed microstructure (Appendix B) Step binding mortar layers was predicted using Dvorkin's model, see Fig. 7a. This skeleton framework was further equivalent as a twophase composite with air voids embedded into a continuous matrix. This two-phase composite's stiffness was assumed to be the same as the skeleton framework. The air voids' volume fraction in the two-phase composite was equal to the sum of the volume fractions of the air voids and the coating mortar layers in the PA mix. Based on the stiffness of the composite and the volume fraction of each phase, the matrix's stiffness can be back-calculated. At last, the coating mortar layers were added into the skeleton framework to form the PA mix. This was equivalent to the process that the same volume of mortar was embedded into the matrix of the two-phase composite. As a result, a three-phase composite, consisting of the matrix, coating mortar layers and air voids, was formed, see Fig. 7b. The volume fractions of the coating mortar layers and the air voids in this three-phase composite were the same as those in the PA mix. Based on the stiffness and the volume fraction of each phase, the stiffness of the three-phase composite was calculated. The stiffness of the PA mix, which was equal to the stiffness of the three-phase composite, became known as well.
Since the matrix in the two-phase composite and the threephase composite was continuous, the back-calculation of the matrix's stiffness and the calculation of the stiffness of the threephase composite can be achieved using CBMM. In this study, the SC model was preferred over others due to the following reasons: (a) this model can be easily solved to back-calculate the stiffness of the matrix (comparing to the GSC model and the Differential model) [33]; and (b) it can provide more accurate predictions when the volume fraction of inclusions is lower than 50% (comparing to the Dilute model and the MT model) [34]. Appendix B describes the procedure to back-calculate the matrix's stiffness of a two-phase composite using the SC model.
Up to now, previous sections briefly introduced Dvorkin's model and presented a framework to predict the stiffness of PA mixes. The next part of the research was to evaluate the performance of Dvorkin's model on the basis of laboratory investigation. Using a Dynamic Shear Rheometer (DSR) device, frequency sweep tests were conducted on laboratory prepared mortar specimens to obtain input parameter G b ; while the PA mix specimens' complex Young's modulus E mix , which was used to evaluate the performance of the model, was measured from uniaxial cyclic compressive tests. In this study, in order to evaluate the performance of Dvorkin's model for PA mixes with different skeleton frameworks, two types of PA mix materials (PA mix-A and PA mix-B) were considered. PA mixes were designed according to the Dutch standard [35,36]. For both mixes, the gradation of the aggregates and their respective densities are shown in Table 1, and the content of asphalt binder with a penetration of 70-100 was 4.3%. The density of the binder was assumed as 1030 kg/m 3 .

Materials and tests
Three replicates of PA mix specimens for uniaxial compressive tests were prepared following the AASHTO T 342-11 standard method [37]. A gyratory compactor was firstly used to make PA mix specimens with a size of U150 mm Â H170 mm. In order to generate different skeleton frameworks, different compaction efforts were performed. The number of the gyrations for making PA mix-A specimens was 5 cycles while for PA mix-B, 40 cycles were performed. As a result, a more air voids content of 18% was achieved for PA mix-A whereas the air voids content of PA mix-B was 13.2%. Since the final size of PA mix specimens for stiffness measurement is 100 mm in diameter and 150 mm in height, specimens from the gyratory compactor were further cored and cut.

Mortar specimens
All the volumes of the asphalt binder, filler and sand (<0.5 mm) in the mix were used to make mortar specimens. By normalizing the proportion of each size of the sand with respect to the maximum size, the gradation of the sand in the mortar was obtained, see Table 2. The asphalt binder's content in the mortar was computed as 23% according to the content of asphalt binder (4.3%) and the proportion of the sand ((1-4.3%)Â14% = 13.4%) in the mix.
The dimension of the mortar specimens is shown in Fig. 8a. The middle part, with a diameter of 6 mm and a height of 12 mm, is used to measure the material's properties; while the end parts, where steel rings with a height of 4 mm and a thickness of 1 mm were fixed, are used to clamp the specimens on the testing machine. A special mould was used to make these specimens, see Fig. 8b. It is highlighted that due to the high asphalt binder content, when the mortar material was heated up to around 160°C, it could smoothly flow. Therefore, there was no compaction work applied during the preparation. More details on the preparation of mortar specimens can be found in reference [38].

Frequency sweep test
Frequency sweep tests with a frequency range from 20 Hz to 0.1 Hz were conducted at five different temperatures of À10°C, 4°C, 21°C, 37°C and 54°C. At each temperature, in order to guarantee the mortar behaves as a linear viscoelastic material, a constant small strain from 10le to 200le was used.

Uniaxial cyclic compressive test
Universal Testing Machine (UTM), see Fig. 9, was used to conduct uniaxial cyclic compressive tests at seven different temperatures, À10°C, 4°C, 21°C, 37°C, 45°C, 54°C and 60°C, and six different frequencies of 20 Hz, 10 Hz, 5 Hz, 1 Hz, 0.5 Hz and 0.1 Hz. Forces were measured via the load cell at the top of the specimen and the displacement was measured via three Linear Variable Differential Transformer sensors. During all the tests, a constant strain with an amplitude of 10le was controlled.

Specified values of the input parameters in Dvorkin's model
Using the proposed method in Section 3, the values of the parameters (i.e. mechanical parameters, volumetric and geometric parameters) in the proposed microstructure were determined on the basis of the materials' properties. In this section, these determined values will be presented.

Mechanical parameters
The master curves of the dynamic shear modulus |G b | and the phase angle d b of the mortar at a reference temperature of 21°C were constructed according to the time-temperature superposition principle, see Fig. 10. It can be observed that with the decrease of frequencies, the value of |G b | keeps decreasing to be close to zero, and meanwhile the value of d b increases to be approximately 90°. This can be related to the fact that the volume fraction of the fine particles in mortar is not high enough to form a particle-on-particle skeleton, and thus its behavior is governed by asphalt binder.
The values of other mechanical parameters m b , G p and m p were obtained from the literature [16,39]  to be sensitive to test frequency and temperature. However, due to the difficulty of measuring the radial deformation of a mortar specimen, it was difficult to directly measure m b from laboratory tests. Therefore, a constant value of m b between 0.35 and 0.5 is generally assumed in the literature [16].

Volumetric and geometrical parameters
As mentioned earlier, the values of / v for PA mix-A and PA mix-B were 18% and 13.2%, respectively. According to the content and the density of asphalt binder, and the aggregate gradation and their respective densities, the values of / b and / p for PA mix-A were calculated as 20% and 62%, respectively, while for PA mix-B, these values were 20.8% and 66%, respectively.
Since PA mix-A and PA mix-B contain the same aggregate gradations and the same asphalt binder contents, the values of  Table 3.
Comparing the values of the geometric parameters between PA mix-A and PA mix-B, it can be found that different compaction effort affects individual phases in terms of not only their volumetric properties in the macroscale but also their geometric properties in the microscale. A larger amount of compaction effort decreases the average distance between adjacent particles (or the average thickness of the binding mortar layers). When two particles become closer, a higher proportion of the coating mortar layers overlap each other to form the binding mortar layers. Therefore, with the decrease of the distance between adjacent particles, the radius, as well as the total volume of the binding mortar layers, increases. On the other hand, more compaction effort decreases the average coordination number. This can be explained by the fact that when the particles become closer, less space around one particle is provided for other particles to surround it.

Predicted results of PA mixes' modulus
In Fig. 11, both the experimental results and the predicted values of E mix (represented in terms of dynamic Young's modulus | E mix | and phase angle d mix ) are presented. It can be seen that the experimental values of |E mix | show asymptotic behaviour both at very high frequencies and very low frequencies. In addition, with the decrease of frequencies, the experimental results of d mix   increase initially and decrease after a peak point. It is highlighted here that previous research work [38], concluded that these behaviours of PA mixes are strongly related to the characteristics of the stone-on-stone skeleton in the mix. Fig. 11 also shows a comparison between the experimental and the predicted values. From the plots, it can be seen that Dvorkin's model can predict E mix and d mix quite well for frequencies higher than around 0.1 Hz. This shows that Dvorkin's model works well within specified frequencies.
In order to check the viability of Dvorkin's model, predicted results were compared against commonly used CBMM in Fig. 11. The Differential model was chosen over other models because previous research studies showed that the Differential model provided better predictions for PA mixes [17]. As can be seen in Fig. 11, at higher frequencies, for a loosely compacted PA mix (i.e. PA mix-A), the predicted results of |E mix | and d mix do not show significant differences between Dvorkin's models and the Differential model. Whereas, for a densely compacted PA mix (PA mix-B), Dvorkin's model provides better predictions. This can be related to the fact that Dvorkin's model is developed mainly for densely packed granular materials [30]. Apart from the volumetric properties of each phase, Dvorkin's model also takes into account the stiffening effects of the geometric characteristics of individual particles and mortar layers. Therefore, in comparison to the Differential model, Dvorkin's model can account for the inter-particle interactions more accurately.
The limitation of Dvorkin's model can also be clearly seen from the plots of Fig. 11 that at lower frequencies, the predicted values of |E mix | are lower than the expected values. Moreover, the predicted results of d mix do not show the expected decreasing trend. One of the possible reasons could be that the values of some input parameters, i.e. h 0 , n, a and m b , should change with frequencies/temperatures rather than being constants over the whole frequency range. With this realization, in order to further understand the limitation of Dvorkin's model, the sensitivity analysis of the predicted stiffness on these parameters was conducted. In the following section, the results of the sensitivity analysis are presented.

Sensitivity analysis of predicted results to input parameters
In the previous analysis, the geometrical characteristics of the proposed PA mix's microstructure were considered to remain unchanged. However, this assumption may not be true because the external loads, i.e. the gravity and the applied load, may change the values of geometrical parameters (i. e. h 0 , n and a). With the increase of temperatures and the decrease of frequencies, the mortar layers become softer and easier to deform. Consequently, the distance between adjacent particles, i.e. the value of h 0 , is expected to decrease. Meanwhile, the change of h 0 also induces the change of n and a, which can be easily derived from Eqs. (B.9)-(B.12) and Eq. (B.16).
The above discussion indicates that the use of constant values of h 0 , n and a in the whole frequency range can be one reason for the poor performance of Dvorkin's model in Fig. 11. Additionally, a constant value of m b was assumed in the predictions. As mentioned earlier, the value of m b is frequency-and temperature-dependent. 6.2.1. The sensitivity of predicted E mix on h 0 As discussed above, the value of h 0 is expected to decrease with the decrease of frequencies, therefore, three lower values of h 0 (<0.22 mm) were used to predict E mix in the sensitivity analysis. Table 4 shows the calculated results of n and a for different values of h 0 using Eqs. (B.9)-(B.12) and Eq. (B.16), respectively. It can be seen that with the decrease of h 0 , the values of n and / b_c decrease while the values of a and / b_b increase. The reasons for the changes of n, a, / b_b and / b_c with the value of h 0 have been explained in the previous section (see Section 5.2). Fig. 12 presents the predicted |E mix |-f and d mix -f curves using the values of h 0 , n and a in Table 4. It can be observed that with the decrease of h 0 , PA mixes become stiffer, which is reflected by an increase of predicted |E mix | and a decrease of predicted d mix . This can be explained by the fact that with the decrease of two particles' distance, the total volume of the binding mortar layers, which make the main contribution to the load-bearing capacity of the mix, increases. Meanwhile, when two particles become closer, their interactions become stronger. As can be derived from Eqs. (4) and (5), the stresses at the mortar-particle interface increase with the decrease of h 0 . When the interactions between adjacent particles become stronger, their stiffening effect on the mix becomes more significant.
Further investigations of the plots in Fig. 12 show that the shapes of the predicted |E mix |-f and d mix -f curves do not significantly change with different values of h 0 . At very low frequencies, even when h 0 is close to 0, the predicted values of |E mix | are still much lower than the experimental results. Moreover, the predicted values of d mix do not show a significant decrease. These observations indicate that the predictions are not likely to be accurate even if the values of the geometric parameters change.
Therefore, it can be concluded that the inaccurate assumption of geometric parameters is not the main reason to explain the poor performance of the model at high temperatures/low frequencies.
With this realization, further investigation was conducted to evaluate the influence of another assumed parameter m b on the predicted results.

The sensitivity of predicted E mix on v b
The predicted results of E mix using different values of v b are shown in Fig. 13 1-2m b )) when the value of m b ranges from 0 to 0.5. It can be observed that when m b is smaller than 0.45, the value of S n_b /2G b is not higher than 6. This indicates that when m b is smaller than 0.45, the value of S n_b is at most 12 times of G b . Since G b is much lower than the stiffness of aggregates, the value of S n_b should be much lower than the stiffness of aggregates as well. In this case, the normal stiffness of the two bonded particles (F n /d n ) and the predicted stiffness of the mix are governed by the softer binder layer. Therefore, in Fig. 13, it was observed that when m b is equal to 0.3 and 0.4, the predicted results of |E mix | and d mix follow the same shapes as |G b | and d b (see Fig. 10), respectively. On the contrary, when m b is higher than 0.45, the value of S n_b increases significantly. When m b approaches 0.5, S n_b becomes infinite, which indicates that the binder layer behaves as a rigid body. In this case, when F n is applied, only aggregates deform and thus the predicted modulus of the mix relies on the properties of the aggregates. Since the aggregates are stiffer and behave as elastic materials, the predicted properties of the mix (see Fig. 13) tend to be stiff (high value of |E mix |) and elastic (low values of d mix ) as well.
From the above discussions, it can be concluded that the effect of m b on predicted |E mix | and d mix is related to the assumption that with the increase of m b , the mortar layer behaves like a rigid body.
In this case, when two aggregate particles approach each other, the mortar layer does not have any deformation, while only aggregates deform. However, in a real PA mix, when m b approaches 0.5 at high temperatures or low frequencies, it is expected that the mortar layer behaves more like a viscous material that easily deforms under loading, while the deformation of the stiff aggregates is much less significant. It is obvious that in these two cases, the properties of the mortar layer are different. Therefore, although the predictions with higher values of t b may match with the experimental results numerically, the physical mechanism behind these predictions does not seem to be realistic.

Limitation of Dvorkin's model and possible explanations
Predicted results in Fig. 11 show that Dvorkin's model performs better in predicting the modulus of PA mixes at higher frequencies than CBMM, especially for mixes with densely packing aggregate particles; however, it still fails to provide accurate predictions at lower frequencies. Following sensitivity analyses further reveal that it is impossible (or unreasonable) to improve the accuracy of the predictions at lower frequencies by varying the values of input parameters.
The limitation of Dvorkin's model can be related to the assumption that in a bonded granular material, a load is always transferred through the mortar layers between adjacent particles. This assumption is clearly reflected by Eqs. (2) and (3), where the total deformation of a two-bonded particle system is the sum of the deformation of the particles and the deformation of the mortar layer. This assumption is suitable at low temperatures/high frequencies because, at these conditions, the mortar layer between two particles is stiff enough to transfer a high level of load from one particle to the other. Even if a direct particle-on-particle contact area forms at the same time, since the total particle-mortar contact area is expected to dominate over the particle-on-particle contact area, the overall stiffness of the mix is governed by the behavior of the mortar layers.
On the other hand, at high temperatures/low frequencies, the mortar layers are too soft to effectively transfer loads among adjacent particles. Comparing to the soft mortar layers, direct particleon-particle contacts are supposed to play a leading role in providing the load transfer capacity for the mix. In this case, the assumption in Dvorkin's model is not applicable anymore, and thus the predictions start to differ from the experimental results.

Conclusions
This study presented a methodology to use Dvorkin's model to predict the stiffness of PA mixes. The predicted results were compared to those from commonly used CBMM and the experimental values. In order to improve the accuracy of the predictions at lower frequencies, the sensitivity of the predicted results on geometrical characteristics of the mix and the Poisson's ratio of the binder layer was analysed. In the end, the limitation of Dvorkin's model was  highlighted. Based on the obtained results, the following conclusions can be drawn: In the proposed microstructure model, a PA mix was simulated as an assembly of identical spherical particles that are covered by a uniform thickness of mortar layers and bonded together by binding mortar layers. Geometric parameters were determined in a way that: the radius of the spherical particles was determined from the aggregate gradation; the thickness of the coating mortar layers was determined from the mortar content; the minimum distance between two adjacent particles and the coordination number were determined from air voids content, and the radius of the binding mortar layers was finally determined on the basis of the values of other parameters. Based on the proposed microstructure model, the stiffness of PA mixes was calculated in three steps. At first, the stiffness of the skeleton framework consisting of aggregate particles and the binding mortar layers was predicted using Dvorkin's model. Then the matrix's stiffness of an equivalent two-phase composite, whose stiffness was identical to that of the skeleton framework, was back-calculated using the Self-consistent model. At last, by adding the coating mortar layers in the matrix, the stiffness of PA mixes was calculated using the Self-consistent model again. Dvorkin's model describes a more realistic microstructure for PA type mixes. This model takes several geometrical characteristics of the aggregate particles and the mortar layers into account. Therefore, the utilization of Dvorkin's model is beneficial to better understand how the characteristics in the microscale affect the mechanical properties of a mix. Furthermore, in a wide range of high frequencies, Dvorkin's model can provide accurate predictions not only for loosely compacted PA mixes but also for PA mixes with densely packed aggregate particles. Dvorkin's models showed inaccurate predictions at lower frequencies: the predicted dynamic moduli were significantly lower than the experimental results, and moreover, the predicted phase angle did not show a decreasing trend. These poor performances cannot be significantly improved by varying the values of the geometric parameters. By varying the Poisson's ratio of the mortar layers, the predictions may match the experimental results numerically. However, the physical mechanism behind these predictions that the mortar layers behave like a rigid body does not seem to be realistic. The limitation of Dvorkin's model is related to the assumption that in a bonded granular material, a load is always transferred through the mortar layers between adjacent particles. This assumption is valid at higher frequencies, while at lower frequencies, the load is supposed to be mainly transferred through the direct particle-on-particle contacts.
In future research work, the capability of Dvorkin's model to accurately predict the stiffness of PA mixes at higher frequencies will be further validated using PA mixes with different gradations and binder contents. Furthermore, Dvorkin's model will be modified to take into account the contribution from the direct particle-on-particle contacts. This modification may be achieved by using Hertzian contact theory to describe the mechanical responses of two contacting particles.

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.

Acknowledgements
The corresponding author thanks the financial support from the China Scholarship Council (CSC).

Data availability
The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.

Appendix A. Numerical solutions to the displacement of binder layer
In this study, a numerical method was used to calculate V(x) and U(x). In the following paragraphs, the steps for solving V(x) in Eq. (9) are presented and the same method can be used to solve U(x) in Eq. (10). It is noted that the method used in this study was developed on the basis of the basic approaches for solving the Volterra integral equations of the second kind [40].
The calculation of V(x) can be made in three steps: Step1: choose the points where the value of V(x) is calculated. In the range from 0 to a, N + 1 discrete points of x with an identical distance between neighboring points, x 0 , x 1 ,. . ., x N , are chosen. By substituting these points into Eq. (9), N + 1 equations are obtained, see Eq. (A.1). The following steps aim to solve these N + 1 equations to obtain the values of V(x) at different points of x.
ds; i ¼ 0; 1; :::; N ðA:1Þ Step 2: discretize the variables. The variable u is identically discretized into N + 1 discrete points in the range from 0 to p, Vð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x 2 i þ s 2 j À 2x i s j cosu k q Þ ¼ Vðx j Þ; i ¼ 0; 1; :::; N ðA:2Þ x 2 i þ s 2 j À 2x i s j cosu k ¼ x 2 j ; j ¼ 0; 1; :::; n i ðA:3Þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi x 2 j À x 2 i sin 2 u k q ; j ¼ 0; 1; :::; n i ðA:4Þ Step while they are equal to 0.5 for the points at the boundaries. h u is the distance between two neighboring points for u, which is equal to p/N. h ij s is the distance between two neighboring points for s, which can be automatically known once the values of s j are obtained. In the proposed microstructure in Fig. 5b, five geometric parameters, i.e. R, t, h 0 , n and a, are required to predict the mix's effective modulus. In this appendix, the proposed methods to determine these parameters are presented. Once the geometric parameters are known, the volume fractions of the binding mortar layers / b_b and the coating mortar layers / b_c can be further determined. Thus, the proposed method to determine / b_b and / b_c is introduced as well.

B.1 Calculation of R
In Dvorkin's model, aggregate particles are modelled as identical spheres. However, in an asphalt mix, the size of the aggregate particles is generally not uniform but graded. Therefore, in order to specify the value of R, it is necessary to obtain a representative size of the aggregate particles in the mix. For this purpose, a commonly used mathematical method of averaging different sizes of graded aggregate particles was utilized [41], see Eqs. (B.1) and (B.2): where P i+1 and P i are the percentages passing the sieve i + 1 and sieve i by the total weight of aggregates, respectively; q i is the density of aggregates retained on sieve i; / i is the volume fraction of aggregate particles retained on sieve i; N is the aggregates' total grades by sieving, i.e., 0.063 mm, 2 mm, 5.6 mm, etc.; and d i+1 and d i are the diameters of sieve i + 1 and sieve i, respectively.

B.2 Calculation of t
The value of t can be determined according to the binder content b in the mix, which is defined as where V b and / b denote the volumes and the volume fraction of the mortar layers, respectively. Assuming that individual particles are surrounded by mortar layers with identical thicknesses, the value of V b can be calculated using Eq. (B.4).
Since the value of b is known from the values of / b and / p , the value of t can be computed from Eq. (B.5).

B.3 Calculation of h 0
As mentioned earlier, the value of h 0 is related to the compaction effort. Since the compaction effort can be reflected by the air voids content, the value of h 0 can be determined according to the air voids content of the mix.
where V total denotes the total volume of the specimen; V mass and W mass denote the volume and the weight of the real mass of the specimen; and q max is the maximum density of the mix, which can be obtained from the density of each constituent. Before any compaction effort is applied (when N = 0), the value of h 0 is expected to be equal to t. In this case, the corresponding value of / v c is 0.22, see Fig. B1. With the increase of the compaction effort, the value of h 0 decreases and thus the value of / v c decreases as well. After conducting a large amount of compaction effort, h 0 is expected to be gradually close to 0, and meanwhile, / v c shows a gradually stable value of 0.13. It is highlighted here that in Fig. B1, there is no point where the tangent to the curve is horizontal; however, by fitting the curve using an exponential function (/ v c ¼ ae ÀbN þ c; where a, b and c are fitting parameters), it was obtained that the minimum value of / v c was 0.13 (c = 0.13). Therefore, it was stated that / v c gradually reached a stable value of 0.13.
Using the values of h 0 and / v c in the above mentioned two extreme cases, the relationship between h 0 /t and / v c was assumed to be expressed as a simple linear function, see Fig. B2. From this relationship, once the air voids content of the PA mix specimen is known, the value of h 0 /t can be obtained.
Two points are worth noting here: For PA mix specimens with different gradations and binder contents, the relationship between h 0 /t and / v c should be different since the packing conditions of the aggregates are different. Therefore, in order to build such a relationship for a certain type of PA mix, there is no need to compact more samples with different gradations. A linear function may not be the most suitable one to describe the relationship between h 0 /t and / v c . However, since h 0 is a geometric parameter in the microscale, it is difficult to determine other values of h 0 during the compaction process. Sophisticated image processing techniques may be a promising method to solve this issue; however, such work is beyond the scope of this study. Furthermore, the agreement between the experimental results and the predictions, which will be shown in the later section, verifies that using such a simple linear function is acceptable.

B.4 Calculation of n
For a PA mix system in Fig. 5b where spherical particles are covered by mortar layers, it is difficult to directly determine its average coordination number n. However, for a ''dry" packing system consisting of only identical spheres, researchers [42] have proposed that its coordination number n d can be estimated on the basis of the air voids content of the system / d , see Fig. B3. The rela-tionship between n d and / d can be developed by plotting their values in different regular packings, i.e. a simple cubic packing (n d = 6, / d = 0.4764), a tetragonal sphenoidal packing (n d = 8, / d = 0.3954), a pyramidal packing (n d = 10, / d = 0.3019) and a tetrahedral packing (n d = 12, / d = 0.2595). As shown in Fig. B3, the value of n d shows an almost linear increase with the decrease of / d . Therefore, a regressed linear equation, see Eq. (B.9), can be given to describe the relationship between n d and / d . / 0 ¼ À0:0374n þ 0:6956 ðB:9Þ On the basis of the relationship in Fig. B3, researchers [42] further developed a method to determine the average coordination number of a packing system where the binding mortar layers were included. Using a similar method, the value of n was determined in this study.
In order to determine the value of n, a reference packing system is introduced, see Fig. B4. In the reference system, individual particles have the same configurations as those in the PA mix in Fig. 5b. This indicates that the coordination number in the reference system remains the same as n. However, particles in the reference system are covered by mortar layers with a thickness of h 0 instead of a thickness of t. Due to the change of the mortar layers' thickness, the air voids content in the reference system changes to a value of / 0 .
Since the minimum distance between adjacent particles in the PA mix system is 2 h 0 , particles covered with a mortar layer of h 0 in the reference system are supposed to just touch each other. Therefore, the arrangement of the mortar-coated particles in the reference system can be considered as the same as a dry packing system consisting of spheres with a uniform radius of (R + h 0 ). Accordingly, the relationship between n and / 0 can be described using the relationship in Eq. (B.9).
Total volumes of the particles and the mortar layers V a ' in the reference system (see Fig. B4) can be calculated using Eq. (B.10), and thus the total volume of the reference system can be given as V a ' /(1-/ 0 ). It is further assumed that more mortar materials are located around the particles in Fig. B4 to form the PA mix system in Fig. 5b. Since the total volume of the system cannot change, the relationship in Eq. (B.11) must satisfy.  where / v is the air voids content of the PA mix.
/ r ðK r À K m Þð3K eff þ 4G eff Þ 3K r þ 4G eff ðC:1Þ 5/ r G eff ðG r À G m Þð3K eff þ 4G eff Þ 3K eff ð3G eff þ 2G r Þ þ 4G eff ð2G eff þ 3G r Þ ðC:2Þ where and the subscripts ''m", ''r" and ''eff" denote the matrix phase, the inclusion phase r and the effective medium, respectively;/ denotes each phase's volume fraction. For the calculation of the stiffness of the equivalent two-phase composite in Fig. 7a where / 2 is equal to the addition of / v and / b c ; and K eff and G eff are identical to the moduli of the skeleton framework in Fig. 7a. Once the values of K eff and G eff are known, the values of K m and G m can be easily derived as: Furthermore, in order to calculate the stiffness of PA mixes as a three-phase composite in Fig. 7b, Eqs. (C.1) and (C.2) can be rearranged as: where / 2 and / 3 are equal to / b c and / v , respectively; and K 2 and G 2 are equal to K b and G b , respectively. It is noted here that it is difficult to obtain explicit solutions for K eff and G eff from Eqs. (C.7) and (C.8); thus, a numerical method, i.e. the Newton-Raphson method, was applied in this study. The details about how to use this method can be found elsewhere [38].