A Molecular Dynamics Approach to the Impacts of Oxidative Aging on the Engineering Characteristics of Asphalt

Oxidative aging is an inevitable environmental factor that accelerates asphalt pavement deterioration. This study employed a molecular dynamics simulation to investigate the impact of aging on asphalt cement from the perspectives of thermodynamic properties, and diffusion and adhesion characteristics. Results indicate that aging increased bulk density from 1.008 to 1.081 g/cm3 and cohesive energy density by 15.6%, which was attributed to the promoted molecular polarity and intermolecular attractiveness. The enhanced molecular interactions also reduced molecular mobility, which led to an increase in the glass transition temperature by 30 K, suggesting that aging diminished the resistance of asphalt to thermal cracking. Simulations of the diffusion behaviors across different temperatures demonstrated that the Arrhenius relationship described well the temperature dependence of the diffusion coefficient, and that aging considerably slowed down the diffusion process as represented by Arrhenius prefactor D0, which dropped by 38.2%. The asphalt–aggregate adhesion was assessed using layered models with and without a water interlayer of different thicknesses. The adhesion was enhanced upon aging due to the significantly improved electrostatic interactions at the interface. Evaluation of the residual adhesion with the presence of interfacial water suggested that aging would raise the moisture susceptibility of asphalt pavement. The increase in molecular polarity was considered to be highly responsible for these aging consequences, and was thus further investigated via the electrostatic potential surface and dipole moment.


Introduction
As a byproduct of crude-oil distillation, asphalt has been most widely used in pavement engineering. It plays a central role in maintaining the material integrity and strength of pavement materials against traffic loading and environmental weathering. Despite ubiquity and macroscopic homogeneity, the chemical composition of asphalt still remains a mystery, as it comprises numerous molecules of different elements, structures, and properties. These molecules are usually categorized and separated as per solubility into the four fractions of saturates, aromatics, resins, and asphaltenes (i.e., the so-called SARA fractionation). There is a continuing debate on the representation of the internal structure of asphalt. Two philosophies primarily exist: one considers asphalt to have a colloidal structure, whereas the other views it as a homogeneous liquid without colloidal characteristics [1]. Nevertheless, to a great extent, asphalt can be practically treated as a mixture of polymers, as many concepts (e.g., glass transition temperature, molecular weight distribution, and complex modulus) and modeling theories (such as viscoelasticity and free volume theory) that are largely rooted in polymer materials have been widely applied to characterize asphalt, and are used in the inference of engineering performance [2][3][4].
As an organic binding agent, asphalt cement is prone to oxidative aging, which leads to degradation in engineering properties such as stress relaxation and healing, thereby

Material Models and Simulation Method
The MD simulation was performed using commercial software Materials Studio [23]. The Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies (COMPASS) II forcefield was employed to describe the molecular interactions. The parameters in this forcefield were derived through a combination of fits to the ab initio data for valence parameters, and experimental data for nonbond parameters [24]. For vdW interactions, the atom-based summation method was used with a cut-off distance of 15.5 Å, and the long-range correction was applied beyond this range. For electrostatic interactions, the Ewald summation was adopted. All the computations were performed with a time step of 1.0 fs.

Asphalt Bulk Models
Asphalt is an extremely complex mixture of hydrocarbons for which the simulation is usually carried out using a set of molecular models representative of the different fractions. Given considerations on the computational cost, the virgin (unaged) asphalt system in this study was established using one molecular structure for each of the four fractions on the basis of Li and Greenfield's model [11]. Specifically, the virgin system consisted of asphaltene-pyrrole, quinolinohopane, perhydrophenanthrene-naphthalene (PHPN), and hopane at a number ratio of 2:6:8:2, as shown in Figure 1a and Table 1. 15.5 Å, and the long-range correction was applied beyond this range. For electrostatic interactions, the Ewald summation was adopted. All the computations were performed with a time step of 1.0 fs.

Asphalt Bulk Models
Asphalt is an extremely complex mixture of hydrocarbons for which the simulation is usually carried out using a set of molecular models representative of the different fractions. Given considerations on the computational cost, the virgin (unaged) asphalt system in this study was established using one molecular structure for each of the four fractions on the basis of Li and Greenfield's model [11]. Specifically, the virgin system consisted of asphaltene-pyrrole, quinolinohopane, perhydrophenanthrene-naphthalene (PHPN), and hopane at a number ratio of 2:6:8:2, as shown in Figure 1a and Table 1.
As previously mentioned, ketones and sulfoxides are the major oxidation products occurring at the benzyl carbons and alkyl/aryl sulfides, respectively. Therefore, the structures for the virgin asphaltene, resin, and aromatic fractions were modified accordingly to represent the oxidized counterparts (Figure 1b) [17]. The nonpolar saturate fraction rarely changes with time [25]; thus, the same structure was retained. In the aged model, the fraction ratio was slightly adjusted to account for the component shift typically observed during the asphalt oxidation process [5,26]. Table 1 provides a compositional summary for the virgin and aged asphalt models.    As previously mentioned, ketones and sulfoxides are the major oxidation products occurring at the benzyl carbons and alkyl/aryl sulfides, respectively. Therefore, the structures for the virgin asphaltene, resin, and aromatic fractions were modified accordingly to represent the oxidized counterparts ( Figure 1b) [17]. The nonpolar saturate fraction rarely changes with time [25]; thus, the same structure was retained. In the aged model, the fraction ratio was slightly adjusted to account for the component shift typically observed during the asphalt oxidation process [5,26]. Table 1 provides a compositional summary for the virgin and aged asphalt models.

Mineral Substrate
Calcite (commonly found in limestone) was used as the mineral substate to evaluate the adhesion properties of the virgin and aged asphalts. The stone aggregates used in asphalt mixtures are obtained mainly through a mechanical crushing process during which the mineral tends to rupture along the surfaces where interlayer bonding strengths are weak. For this reason, the crystal first needs to be cleaved to expose such surfaces using the Morphology module in Materials Studio. Figure 2a depicts the morphology calculation result based on attachment energy theory [27]. This theory considers attachment energy as a measure of crystal growth rate normal to a face, and that the face with higher attachment energy would grow faster and thus be less morphologically important. The use of this theory resulted in a single family of lattice planes that were represented in terms of Miller indices as (1 0 4). This family of planes was also identified using the Bravais-Friedel-Donnay-Harker (BFDH) approach [28], but with the second highest morphological importance (not shown). BFDH theory only considers the lattice geometry in computation without taking into account the details such as chemical nature and molecular packing in the crystal.

Mineral Substrate
Calcite (commonly found in limestone) was used as the mineral substate to evaluate the adhesion properties of the virgin and aged asphalts. The stone aggregates used in asphalt mixtures are obtained mainly through a mechanical crushing process during which the mineral tends to rupture along the surfaces where interlayer bonding strengths are weak. For this reason, the crystal first needs to be cleaved to expose such surfaces using the Morphology module in Materials Studio. Figure 2a depicts the morphology calculation result based on attachment energy theory [27]. This theory considers attachment energy as a measure of crystal growth rate normal to a face, and that the face with higher attachment energy would grow faster and thus be less morphologically important. The use of this theory resulted in a single family of lattice planes that were represented in terms of Miller indices as (1 0 4). This family of planes was also identified using the Bravais-Friedel-Donnay-Harker (BFDH) approach [28], but with the second highest morphological importance (not shown). BFDH theory only considers the lattice geometry in computation without taking into account the details such as chemical nature and molecular packing in the crystal.
Considering additional experimental evidence confirming that surface (1 0 4) is the most commonly exposed plane for calcite [29,30], this habit face was used for crystal cleavage in this study. Supercells with appropriate dimensions depending on the need (discussed later) were then constructed for the interface models. A fractional thickness of 6 was applied, and the resulting supercell substrates had a thickness of 17 Å (greater than the vdW cut-off distance), as illustrated in Figure 2b.

General Setup for the MD Simulation
For constructing the bulk model of asphalt, the Amorphous Cell (AC) module was utilized, which allowed for creating cells with periodical boundary conditions containing random arrangements of the model molecules. The initial density was set at 0.6 g/cm 3 (lower than the equilibrium level) to encourage the randomness in generating the asphalt AC configurations.
In this study, unless otherwise stated, the completion of the model setup was followed by a four-step process, as illustrated in Figure 3, to reach the final state at equilibrium that is ready for analysis.

•
Step 1: geometry optimization to eliminate the unrealistic molecular overlapping and adjust unstable high-energy configurations. Considering additional experimental evidence confirming that surface (1 0 4) is the most commonly exposed plane for calcite [29,30], this habit face was used for crystal cleavage in this study. Supercells with appropriate dimensions depending on the need (discussed later) were then constructed for the interface models. A fractional thickness of 6 was applied, and the resulting supercell substrates had a thickness of 17 Å (greater than the vdW cut-off distance), as illustrated in Figure 2b.

General Setup for the MD Simulation
For constructing the bulk model of asphalt, the Amorphous Cell (AC) module was utilized, which allowed for creating cells with periodical boundary conditions containing random arrangements of the model molecules. The initial density was set at 0.6 g/cm 3 (lower than the equilibrium level) to encourage the randomness in generating the asphalt AC configurations.
In this study, unless otherwise stated, the completion of the model setup was followed by a four-step process, as illustrated in Figure 3, to reach the final state at equilibrium that is ready for analysis. Step 3: equilibration to bring the system to equilibration over a simulation time of 500 ps using the NVT ensemble with an Andersen thermostat.

•
Step 4: data production run to further equilibrate the system under the NVT or NPT ensemble (a constant number of atoms N, a controlled pressure P, and a controlled temperature T) depending on the model of interest, using the Nosé-Hoover thermostat (for physically sound dynamics) and Berendsen barostat, over a simulation time of 1 ns. Figure 3. Four-step process employed to bring the system to equilibrium with physically sound dynamics (time not to scale).

Bulk Thermodynamic Properties
This section presents the thermodynamic properties of the virgin and aged bulk asphalts predicted from the MD simulation, including density, cohesive energy density (CED), and glass transition temperature Tg. Results preliminarily validated the established molecular model and the MD simulation process by comparing the existing experimental and simulation findings.

Prediction of Density
The four-component asphalt models were brought to equilibrium under atmospheric pressure (1 atm) and room temperature (298.15 K) using the four-step simulation process (NPT for the last step). The final density for the virgin asphalt was equilibrated at 1.008 g/cm 3 , which agreed well with the existing simulation results and was comparable with experimental values in the range from 1.0 to 1.04 g/cm 3 [11,17,20]. The aged model yielded a slightly higher density of 1.081 g/cm 3 . The increase in density could be partly attributed to the fact that oxidation increased the molecular weight of the polar fractions by introducing oxygen atoms and tilted the composition balance towards asphaltenes, which was reflected in the molecular structures and proportions used in simulation.

Prediction of Cohesive Energy Density
The CED of a material in a condensed state is defined as the amount of energy required to completely remove a unit volume of molecules against all attractive forces from the neighboring molecules. Hence, it can be used as an overall measure of the strength of intermolecular interactions. The predicted CED values for the virgin and aged models (eventually equilibrated under 1 atm) were 3.02 × 10 8 and 3.49 × 10 8 J/m 3 , respectively. The primary contribution to CED was vdW interactions, and aging slightly reduced the vdW contribution from 96.7% to 91.9% due to the increased molecular polarity. These observations agree with values reported in the literature [15,31].

•
Step 1: geometry optimization to eliminate the unrealistic molecular overlapping and adjust unstable high-energy configurations.

•
Step 2: annealing to overcome energy barriers by periodically heating and cooling the system (250 to 550 K) to find the energetically favorable minima using the canonical ensemble NVT (a constant number of atoms N, a constant volume V, and a controlled temperature T).

•
Step 3: equilibration to bring the system to equilibration over a simulation time of 500 ps using the NVT ensemble with an Andersen thermostat.

•
Step 4: data production run to further equilibrate the system under the NVT or NPT ensemble (a constant number of atoms N, a controlled pressure P, and a controlled temperature T) depending on the model of interest, using the Nosé-Hoover thermostat (for physically sound dynamics) and Berendsen barostat, over a simulation time of 1 ns.

Bulk Thermodynamic Properties
This section presents the thermodynamic properties of the virgin and aged bulk asphalts predicted from the MD simulation, including density, cohesive energy density (CED), and glass transition temperature T g . Results preliminarily validated the established molecular model and the MD simulation process by comparing the existing experimental and simulation findings.

Prediction of Density
The four-component asphalt models were brought to equilibrium under atmospheric pressure (1 atm) and room temperature (298.15 K) using the four-step simulation process (NPT for the last step). The final density for the virgin asphalt was equilibrated at 1.008 g/cm 3 , which agreed well with the existing simulation results and was comparable with experimental values in the range from 1.0 to 1.04 g/cm 3 [11,17,20]. The aged model yielded a slightly higher density of 1.081 g/cm 3 . The increase in density could be partly attributed to the fact that oxidation increased the molecular weight of the polar fractions by introducing oxygen atoms and tilted the composition balance towards asphaltenes, which was reflected in the molecular structures and proportions used in simulation.

Prediction of Cohesive Energy Density
The CED of a material in a condensed state is defined as the amount of energy required to completely remove a unit volume of molecules against all attractive forces from the neighboring molecules. Hence, it can be used as an overall measure of the strength of intermolecular interactions. The predicted CED values for the virgin and aged models (eventually equilibrated under 1 atm) were 3.02 × 10 8 and 3.49 × 10 8 J/m 3 , respectively. The primary contribution to CED was vdW interactions, and aging slightly reduced the vdW contribution from 96.7% to 91.9% due to the increased molecular polarity. These observations agree with values reported in the literature [15,31].

Prediction of Glass Transition Temperature
Asphalt exhibits a glassy state when the temperature is below T g , and rubbery behavior above T g . The glass transition temperature can be thermodynamically interpreted as the threshold below which the reduction in free volume becomes significantly slow over time as temperature decreases [32]. The physical hardening phenomenon (stiffness increase over time) of asphalt occurs at temperatures around T g [33]; more importantly, this parameter is closely related to the propensity of asphalt pavement for thermal cracking [34].
The glass transition temperature is commonly determined as the temperature at which the two asymptotes to the glassy and rubbery regions intersect in the space of specific volume versus temperature. During the simulation process, a trajectory file was first created that contained three different AC configurations as replicates for the purpose of minimizing the potential scattering in T g due to AC randomness. A Perl script was prepared that iterated over every configuration, executing the following algorithm: a short NVT annealing followed by NPT equilibration runs at decreasing temperatures from 600 to 100 K with a step size of 2 K. At each temperature step, the density was checked every 10 ps, and the equilibration dynamics run was performed until the density had stabilized within a tolerance of 0.02 g/cm 3 . Then, a further NPT production run over 10 ps was used to obtain the final density data before entering the next temperature step. Figure 4 presents the simulation data of averaged specific volume versus temperature and T g determination for the virgin and aged asphalts. The virgin system yielded a glass transition temperature of 315 K, which was in line with the typical values reported from existing simulations [20,35]. Compared to the experimental range of 223 to 303 K [32], this T g result was slightly higher, which could be partly attributed to the higher cooling rate used in simulation [36]. The aged asphalt model provided an increased T g value of 345 K, as confirmed by earlier experimental findings [37]. In addition, by comparing the change in the slope of the two asymptotes, the reduction in the free volume of the aged asphalt exhibited an overall lower sensitivity to temperature drop. This observation can be ascribed to the stronger intermolecular attractiveness that restricted the molecular mobility in the aged system.

Prediction of Glass Transition Temperature
Asphalt exhibits a glassy state when the temperature is below Tg, and rubbery behavior above Tg. The glass transition temperature can be thermodynamically interpreted as the threshold below which the reduction in free volume becomes significantly slow over time as temperature decreases [32]. The physical hardening phenomenon (stiffness increase over time) of asphalt occurs at temperatures around Tg [33]; more importantly, this parameter is closely related to the propensity of asphalt pavement for thermal cracking [34].
The glass transition temperature is commonly determined as the temperature at which the two asymptotes to the glassy and rubbery regions intersect in the space of specific volume versus temperature. During the simulation process, a trajectory file was first created that contained three different AC configurations as replicates for the purpose of minimizing the potential scattering in Tg due to AC randomness. A Perl script was prepared that iterated over every configuration, executing the following algorithm: a short NVT annealing followed by NPT equilibration runs at decreasing temperatures from 600 to 100 K with a step size of 2 K. At each temperature step, the density was checked every 10 ps, and the equilibration dynamics run was performed until the density had stabilized within a tolerance of 0.02 g/cm 3 . Then, a further NPT production run over 10 ps was used to obtain the final density data before entering the next temperature step. Figure 4 presents the simulation data of averaged specific volume versus temperature and Tg determination for the virgin and aged asphalts. The virgin system yielded a glass transition temperature of 315 K, which was in line with the typical values reported from existing simulations [20,35]. Compared to the experimental range of 223 to 303 K [32], this Tg result was slightly higher, which could be partly attributed to the higher cooling rate used in simulation [36]. The aged asphalt model provided an increased Tg value of 345 K, as confirmed by earlier experimental findings [37]. In addition, by comparing the change in the slope of the two asymptotes, the reduction in the free volume of the aged asphalt exhibited an overall lower sensitivity to temperature drop. This observation can be ascribed to the stronger intermolecular attractiveness that restricted the molecular mobility in the aged system.

Asphalt Diffusion Characteristics
The potential of asphalt diffusion was evaluated in a setup that simulated the crack healing process [38]. As cracking in reality occurs well after the system reaches thermodynamic equilibrium, the equilibrated AC configuration was used to build the diffusion model as opposed to the use of an unequilibrated system in the literature [39].

Asphalt Diffusion Characteristics
The potential of asphalt diffusion was evaluated in a setup that simulated the crack healing process [38]. As cracking in reality occurs well after the system reaches thermodynamic equilibrium, the equilibrated AC configuration was used to build the diffusion model as opposed to the use of an unequilibrated system in the literature [39]. Specifically, the two models of the virgin or aged asphalt were arranged in a simulation cell with a 5 Å vacuum space in between to represent the crack, as shown in Figure 5. Meanwhile, at the two ends, a 17 Å thick vacuum slab was applied in order to eliminate interference from the mirror images of the model. Afterwards, the MD simulation was performed using the NPT ensemble at 1 atm for a duration of 1 ns to inspect the diffusion process. Specifically, the two models of the virgin or aged asphalt were arranged in a simulation cell with a 5 Å vacuum space in between to represent the crack, as shown in Figure 5. Meanwhile, at the two ends, a 17 Å thick vacuum slab was applied in order to eliminate interference from the mirror images of the model. Afterwards, the MD simulation was performed using the NPT ensemble at 1 atm for a duration of 1 ns to inspect the diffusion process. The diffusion process was measured by the commonly used concept of mean square displacement (MSD), which describes the translational mobility of particles over time. It can be perceived as the spatial extent that is "explored" by the randomly moving molecules due to diffusion. The MSD as a function of time is defined as where ri(t) is the position vector of the i-th particle at time t, ri(0) denotes the reference position at time 0, and N is the total number of particles to be averaged. Figure 6a presents the time histories of MSD for the virgin and aged asphalts at the room temperature of 298.15 K. The virgin system exhibited higher molecular mobility than that of the aged system throughout the whole process. Specifically, at the simulation time of 300 ps, the two models in the virgin system already diffused across the crack space, and the crack was morphologically healed. The two models in the aged system, however, were still in the process of approaching each other with a noticeable void space in between.
(a) (b) Another perspective to examine the diffusion rate relied on the use of interaction energy between the two models in the cell. To account for the difference and variation in the cell dimensions, the interaction energy was further divided by the cross-sectional area of the cell, and results are compared in Figure 6b. The negative energies at the initial states The diffusion process was measured by the commonly used concept of mean square displacement (MSD), which describes the translational mobility of particles over time. It can be perceived as the spatial extent that is "explored" by the randomly moving molecules due to diffusion. The MSD as a function of time is defined as where r i (t) is the position vector of the i-th particle at time t, r i (0) denotes the reference position at time 0, and N is the total number of particles to be averaged. Figure 6a presents the time histories of MSD for the virgin and aged asphalts at the room temperature of 298.15 K. The virgin system exhibited higher molecular mobility than that of the aged system throughout the whole process. Specifically, at the simulation time of 300 ps, the two models in the virgin system already diffused across the crack space, and the crack was morphologically healed. The two models in the aged system, however, were still in the process of approaching each other with a noticeable void space in between. Specifically, the two models of the virgin or aged asphalt were arranged in a simulation cell with a 5 Å vacuum space in between to represent the crack, as shown in Figure 5. Meanwhile, at the two ends, a 17 Å thick vacuum slab was applied in order to eliminate interference from the mirror images of the model. Afterwards, the MD simulation was performed using the NPT ensemble at 1 atm for a duration of 1 ns to inspect the diffusion process. Figure 5. Geometry of the diffusion model.
The diffusion process was measured by the commonly used concept of mean square displacement (MSD), which describes the translational mobility of particles over time. It can be perceived as the spatial extent that is "explored" by the randomly moving molecules due to diffusion. The MSD as a function of time is defined as where ri(t) is the position vector of the i-th particle at time t, ri(0) denotes the reference position at time 0, and N is the total number of particles to be averaged. Figure 6a presents the time histories of MSD for the virgin and aged asphalts at the room temperature of 298.15 K. The virgin system exhibited higher molecular mobility than that of the aged system throughout the whole process. Specifically, at the simulation time of 300 ps, the two models in the virgin system already diffused across the crack space, and the crack was morphologically healed. The two models in the aged system, however, were still in the process of approaching each other with a noticeable void space in between.
(a) (b) Another perspective to examine the diffusion rate relied on the use of interaction energy between the two models in the cell. To account for the difference and variation in the cell dimensions, the interaction energy was further divided by the cross-sectional area of the cell, and results are compared in Figure 6b. The negative energies at the initial states Another perspective to examine the diffusion rate relied on the use of interaction energy between the two models in the cell. To account for the difference and variation in the cell dimensions, the interaction energy was further divided by the cross-sectional area of the cell, and results are compared in Figure 6b. The negative energies at the initial states of both systems suggested that the subsequent diffusion to heal the crack was thermodynamically feasible. A simulation time of 1 ns was sufficient to stabilize the interaction energy for both systems at room temperature, and the aged asphalt reached the energy equilibrium again at a lower rate. At 300 ps, the interaction energy per unit area for the virgin system was around −130 mJ/m 2 , close to the equilibrated level, but the value for the aged system was not much different from the starting point. Existing findings indicated that, for individual molecules, the diffusion rate is dependent on the factors including molecular weight, shape, and polarity [40,41]. Herein, the lower mobility of the aged asphalt as an assembly of oxidized molecules can be attributed to the increase in molecular polarity and weight that promoted the attractive interactions.
The MD simulation on the diffusion model was further performed at a number of higher temperatures: 333.15, 373.15, 423.15, and 493.15 K (i.e., 60, 100, 150, and 220 • C). The diffusion coefficient D for both the virgin and aged asphalts was obtained from the slope of the MSD plot (by a factor of 1/6), and its dependence on temperature was investigated. Figure 7 presents the relationships of D versus temperature and the fitting using the Arrhenius equation [12]: where E a is the activation energy, R is the universal gas constant (8.314 J/mol/K), T is temperature, and D 0 is the diffusion prefactor.
of both systems suggested that the subsequent diffusion to heal the crack was thermodynamically feasible. A simulation time of 1 ns was sufficient to stabilize the interaction energy for both systems at room temperature, and the aged asphalt reached the energy equilibrium again at a lower rate. At 300 ps, the interaction energy per unit area for the virgin system was around −130 mJ/m 2 , close to the equilibrated level, but the value for the aged system was not much different from the starting point. Existing findings indicated that, for individual molecules, the diffusion rate is dependent on the factors including molecular weight, shape, and polarity [40,41]. Herein, the lower mobility of the aged asphalt as an assembly of oxidized molecules can be attributed to the increase in molecular polarity and weight that promoted the attractive interactions. The MD simulation on the diffusion model was further performed at a number of higher temperatures: 333.15, 373.15, 423.15, and 493.15 K (i.e., 60, 100, 150, and 220 °C). The diffusion coefficient D for both the virgin and aged asphalts was obtained from the slope of the MSD plot (by a factor of 1/6), and its dependence on temperature was investigated. Figure 7 presents the relationships of D versus temperature and the fitting using the Arrhenius equation [12]: where Ea is the activation energy, R is the universal gas constant (8.314 J/mol/K), T is temperature, and D0 is the diffusion prefactor.  Figure 7 shows that, at all temperatures, the virgin system diffused at considerably higher rates than the aged system did, and the diffusion coefficients followed the Arrhenius behavior reasonably well. The slope of the fit in the space of logD versus 1/T was used to further obtain the activation energy, which was 4.4 and 4.0 kJ/mol for the virgin and aged asphalts, respectively. The virgin model yielded a higher activation energy for self-healing but the difference was marginal. The prefactor was 1.1 × 10 −5 cm 2 /s for the virgin, considerably higher than the 6.8 × 10 −6 cm 2 /s for the aged. Li and Greenfield [12] evaluated the correlations of activation energy and prefactor with respect to molecular weight. They reported that the activation energy for different asphalt molecules varied within a narrow range and had little dependence on the molecular weight, but the prefactor showed a clear decreasing trend with it. Considering the fact that the molecular weight increases with aging, prefactor D0 appeared to be a valid index representing the overall aging effect on the asphalt diffusion behavior. In addition, since the increased molecular polarity is also a pronounced aging outcome, further effort is worthwhile to examine its correlation with the diffusion parameters.  Figure 7 shows that, at all temperatures, the virgin system diffused at considerably higher rates than the aged system did, and the diffusion coefficients followed the Arrhenius behavior reasonably well. The slope of the fit in the space of logD versus 1/T was used to further obtain the activation energy, which was 4.4 and 4.0 kJ/mol for the virgin and aged asphalts, respectively. The virgin model yielded a higher activation energy for selfhealing but the difference was marginal. The prefactor was 1.1 × 10 −5 cm 2 /s for the virgin, considerably higher than the 6.8 × 10 −6 cm 2 /s for the aged. Li and Greenfield [12] evaluated the correlations of activation energy and prefactor with respect to molecular weight. They reported that the activation energy for different asphalt molecules varied within a narrow range and had little dependence on the molecular weight, but the prefactor showed a clear decreasing trend with it. Considering the fact that the molecular weight increases with aging, prefactor D 0 appeared to be a valid index representing the overall aging effect on the asphalt diffusion behavior. In addition, since the increased molecular polarity is also a pronounced aging outcome, further effort is worthwhile to examine its correlation with the diffusion parameters.

Asphalt-Aggregate Adhesion Characteristics
In this section, the asphalt-aggregate adhesion property is considered in both dry and wet conditions. The potential size effect was first investigated to determine the proper geometry of the mineral substrate used in simulation, for a balance between modeling accuracy and computational efficiency.

Size Effect
Due to the randomness in generating the asphalt AC, different atoms and functional groups (especially the carbonyl in the aged case) appear at the asphalt-aggregate interface, which may affect the determination of adhesion energy. Statistically, this effect should diminish as the model size enlarges. For this reason, the size effect was first evaluated using the aged asphalt model to determine the proper dimension for the calcite substrate as previously shown in Figure 2b.
In this process, the supercell substrate was first prepared with three different geometries: 32 × 30, 48 × 45, and 64 × 60 Å 2 , as indicated in Figure 8. For each geometry, four different asphalt AC configurations were constructed at the already-known equilibrium density by matching the interface geometry while maintaining the thickness as greater than the vdW cut-off distance. The interface model was created by fitting the asphalt AC onto the supercell substrate, while placing a 50 Å thick vacuum slab on top to avoid interference from the mirror images. Each of the 12 resulting interface models was then subjected to the four-step simulation process described earlier in Figure 3. In the substrate, the three bottom layers of atoms were fixed in their Cartesian coordinates, and the NVT ensemble was used in the final data production step.
Due to the randomness in generating the asphalt AC, different atoms and functional groups (especially the carbonyl in the aged case) appear at the asphalt-aggregate interface, which may affect the determination of adhesion energy. Statistically, this effect should diminish as the model size enlarges. For this reason, the size effect was first evaluated using the aged asphalt model to determine the proper dimension for the calcite substrate as previously shown in Figure 2b.
In this process, the supercell substrate was first prepared with three different geometries: 32 × 30, 48 × 45, and 64 × 60 Å 2 , as indicated in Figure 8. For each geometry, four different asphalt AC configurations were constructed at the already-known equilibrium density by matching the interface geometry while maintaining the thickness as greater than the vdW cut-off distance. The interface model was created by fitting the asphalt AC onto the supercell substrate, while placing a 50 Å thick vacuum slab on top to avoid interference from the mirror images. Each of the 12 resulting interface models was then subjected to the four-step simulation process described earlier in Figure 3. In the substrate, the three bottom layers of atoms were fixed in their Cartesian coordinates, and the NVT ensemble was used in the final data production step. The adhesion energy, described in Equation (3), was calculated for each scenario to obtain the mean and coefficient of variation (CoV), as included in Figure 8. The adhesion energy values were comparable to the experimental results reported in the literature [19,42]. The averaged adhesion energy appeared to not have been significantly affected by the three geometries selected for inspection. Nevertheless, variability due to AC randomness presented a decreasing trend according to CoV, as noted in earlier studies [43]. The substrate geometry of 64 × 60 Å 2 yielded the lowest variability and was thus adopted for subsequent calculations.

Dry Condition
The dry condition was evaluated to provide a baseline for the ensuing analysis involving the presence of water. In this case, the interface model consisted of the calcite  (3), was calculated for each scenario to obtain the mean and coefficient of variation (CoV), as included in Figure 8. The adhesion energy values were comparable to the experimental results reported in the literature [19,42]. The averaged adhesion energy appeared to not have been significantly affected by the three geometries selected for inspection. Nevertheless, variability due to AC randomness presented a decreasing trend according to CoV, as noted in earlier studies [43]. The substrate geometry of 64 × 60 Å 2 yielded the lowest variability and was thus adopted for subsequent calculations.

Dry Condition
The dry condition was evaluated to provide a baseline for the ensuing analysis involving the presence of water. In this case, the interface model consisted of the calcite supercell in the geometry of 64 × 60 Å 2 and an AC for the virgin or aged asphalt. The adhesion energy for the asphalt-aggregate interface was determined in a way similar to the concept of interaction energy per unit area used previously in the discussion of diffusion. Specifically, the adhesion energy was calculated as where E asph-agg is adhesion energy, A is interface area, E int is interaction energy, E interf is the total potential energy of the whole interface model, and E asph and E agg are potential energies of the asphalt AC and mineral substrate, respectively. Negative values of E asph-agg indicate attractive interactions. For facilitating data presentation, hereinafter the adhesion energies were given in terms of magnitude by dropping the negative sign, as shown in Figure 9.
Polymers 2022, 14, x FOR PEER REVIEW 10 of 15 supercell in the geometry of 64 × 60 Å 2 and an AC for the virgin or aged asphalt. The adhesion energy for the asphalt-aggregate interface was determined in a way similar to the concept of interaction energy per unit area used previously in the discussion of diffusion. Specifically, the adhesion energy was calculated as where Easph-agg is adhesion energy, A is interface area, Eint is interaction energy, Einterf is the total potential energy of the whole interface model, and Easph and Eagg are potential energies of the asphalt AC and mineral substrate, respectively. Negative values of Easph-agg indicate attractive interactions. For facilitating data presentation, hereinafter the adhesion energies were given in terms of magnitude by dropping the negative sign, as shown in Figure 9.  Figure 9 also includes the vdW and electrostatic components of the total adhesion energies. Oxidation had a pronounced effect in improving the adhesion between asphalt and calcite. In the virgin case, the interface adhesion was primarily attributed to the vdW forces. After aging, the vdW contribution slightly dropped, but the electrostatic component was substantially increased to a level comparable to the vdW part. The raised importance of electrostatic contribution was ascribed to the interactions between the positively charged calcium at the mineral surface and the strong electronegativity of oxygen introduced into the aged asphalt. These observations were in line with general findings from experiments and simulations [15,44].

Wet Condition
The wet condition was considered to be water that had penetrated into the interface and was present as an interlayer. The severity of moisture damage was simulated using four different water layer thicknesses: 3, 6, 9, and 12 Å, consisting of 400, 800, 1200, and 1600 water molecules, respectively. The previously detailed four-step process was applied to equilibrate the obtained three-layer interface models at the temperature of 298.15 K. Figure 10 shows the model states before and after equilibration for both the virgin and aged cases with a 3Å water interlayer.   Figure 9 also includes the vdW and electrostatic components of the total adhesion energies. Oxidation had a pronounced effect in improving the adhesion between asphalt and calcite. In the virgin case, the interface adhesion was primarily attributed to the vdW forces. After aging, the vdW contribution slightly dropped, but the electrostatic component was substantially increased to a level comparable to the vdW part. The raised importance of electrostatic contribution was ascribed to the interactions between the positively charged calcium at the mineral surface and the strong electronegativity of oxygen introduced into the aged asphalt. These observations were in line with general findings from experiments and simulations [15,44].

Wet Condition
The wet condition was considered to be water that had penetrated into the interface and was present as an interlayer. The severity of moisture damage was simulated using four different water layer thicknesses: 3, 6, 9, and 12 Å, consisting of 400, 800, 1200, and 1600 water molecules, respectively. The previously detailed four-step process was applied to equilibrate the obtained three-layer interface models at the temperature of 298.15 K. Figure 10 shows the model states before and after equilibration for both the virgin and aged cases with a 3Å water interlayer. supercell in the geometry of 64 × 60 Å 2 and an AC for the virgin or aged asphalt. The adhesion energy for the asphalt-aggregate interface was determined in a way similar to the concept of interaction energy per unit area used previously in the discussion of diffusion. Specifically, the adhesion energy was calculated as where Easph-agg is adhesion energy, A is interface area, Eint is interaction energy, Einterf is the total potential energy of the whole interface model, and Easph and Eagg are potential energies of the asphalt AC and mineral substrate, respectively. Negative values of Easph-agg indicate attractive interactions. For facilitating data presentation, hereinafter the adhesion energies were given in terms of magnitude by dropping the negative sign, as shown in Figure 9.  Figure 9 also includes the vdW and electrostatic components of the total adhesion energies. Oxidation had a pronounced effect in improving the adhesion between asphalt and calcite. In the virgin case, the interface adhesion was primarily attributed to the vdW forces. After aging, the vdW contribution slightly dropped, but the electrostatic component was substantially increased to a level comparable to the vdW part. The raised importance of electrostatic contribution was ascribed to the interactions between the positively charged calcium at the mineral surface and the strong electronegativity of oxygen introduced into the aged asphalt. These observations were in line with general findings from experiments and simulations [15,44].

Wet Condition
The wet condition was considered to be water that had penetrated into the interface and was present as an interlayer. The severity of moisture damage was simulated using four different water layer thicknesses: 3, 6, 9, and 12 Å, consisting of 400, 800, 1200, and 1600 water molecules, respectively. The previously detailed four-step process was applied to equilibrate the obtained three-layer interface models at the temperature of 298.15 K. Figure 10 shows the model states before and after equilibration for both the virgin and aged cases with a 3Å water interlayer.  The adhesion energies of interest included the one between water and calcite E water-agg , and the one between water and asphalt E water-asph . In addition, despite the presence of water, there were still some weak interactions existing between asphalt and calcite, and this residual adhesion energy E asph-agg,res was also calculated. Figure 11 shows these energy results as a function of the water layer thickness. Calcite is hydrophilic in nature, and E water-agg was expectedly much higher than the remaining two in Figure 11; it was also an order of magnitude greater than the asphalt-calcite adhesion in the dry case as previously shown in Figure 9. The asphalt-calcite adhesion was significantly reduced with the presence of water, and diminished further with increase in the water layer thickness. Aging consistently improved the interactions of asphalt with water and with calcite at all water layer thicknesses according to E water-asph and E asph-agg,res , respectively. This improvement was attributable to the electronegativity of the added oxygens and thus the increased molecular polarity.
The adhesion energies of interest included the one between water and calcite Ewateragg, and the one between water and asphalt Ewater-asph. In addition, despite the presence of water, there were still some weak interactions existing between asphalt and calcite, and this residual adhesion energy Easph-agg,res was also calculated. Figure 11 shows these energy results as a function of the water layer thickness. Calcite is hydrophilic in nature, and Ewateragg was expectedly much higher than the remaining two in Figure 11; it was also an order of magnitude greater than the asphalt-calcite adhesion in the dry case as previously shown in Figure 9. The asphalt-calcite adhesion was significantly reduced with the presence of water, and diminished further with increase in the water layer thickness. Aging consistently improved the interactions of asphalt with water and with calcite at all water layer thicknesses according to Ewater-asph and Easph-agg,res, respectively. This improvement was attributable to the electronegativity of the added oxygens and thus the increased molecular polarity.
A comparison of the adhesion energies before and after water penetration in terms of the Easph-agg,res/Easph-agg ratio was proposed as an indicator for the moisture susceptibility of the asphalt, as shown in Figure 12. The aged asphalt yielded lower ratios at all the four severity levels. This simulation result led to the inference that oxidation rendered the asphalt more sensitivity to moisture damage, which agrees with the general observations in the asphalt pavement industry.

Further Inspection of Increased Molecular Polarity
A review of the preceding investigations revealed that an increase in the molecular polarity played a critical role in nearly all the aging-induced changes of the evaluated asphalt properties. The elevated polarity was largely responsible for the enhanced A comparison of the adhesion energies before and after water penetration in terms of the E asph-agg,res /E asph-agg ratio was proposed as an indicator for the moisture susceptibility of the asphalt, as shown in Figure 12. The aged asphalt yielded lower ratios at all the four severity levels. This simulation result led to the inference that oxidation rendered the asphalt more sensitivity to moisture damage, which agrees with the general observations in the asphalt pavement industry. The adhesion energies of interest included the one between water and calcite Ewateragg, and the one between water and asphalt Ewater-asph. In addition, despite the presence of water, there were still some weak interactions existing between asphalt and calcite, and this residual adhesion energy Easph-agg,res was also calculated. Figure 11 shows these energy results as a function of the water layer thickness. Calcite is hydrophilic in nature, and Ewateragg was expectedly much higher than the remaining two in Figure 11; it was also an order of magnitude greater than the asphalt-calcite adhesion in the dry case as previously shown in Figure 9. The asphalt-calcite adhesion was significantly reduced with the presence of water, and diminished further with increase in the water layer thickness. Aging consistently improved the interactions of asphalt with water and with calcite at all water layer thicknesses according to Ewater-asph and Easph-agg,res, respectively. This improvement was attributable to the electronegativity of the added oxygens and thus the increased molecular polarity. A comparison of the adhesion energies before and after water penetration in terms of the Easph-agg,res/Easph-agg ratio was proposed as an indicator for the moisture susceptibility of the asphalt, as shown in Figure 12. The aged asphalt yielded lower ratios at all the four severity levels. This simulation result led to the inference that oxidation rendered the asphalt more sensitivity to moisture damage, which agrees with the general observations in the asphalt pavement industry.

Further Inspection of Increased Molecular Polarity
A review of the preceding investigations revealed that an increase in the molecular polarity played a critical role in nearly all the aging-induced changes of the evaluated asphalt properties. The elevated polarity was largely responsible for the enhanced

Further Inspection of Increased Molecular Polarity
A review of the preceding investigations revealed that an increase in the molecular polarity played a critical role in nearly all the aging-induced changes of the evaluated asphalt properties. The elevated polarity was largely responsible for the enhanced intermolecular attractiveness and weakened molecular mobility in the bulk asphalt, and also is the cause for the increased asphalt-aggregate interactions in both the dry and wet conditions. It is, therefore, of great importance and interest to further look into the specific impacts of aging on the molecular polarity. The asphaltene molecules were selected for this purpose since this fraction mostly dominates the engineering properties of asphalt and also presents the highest chemical reactivity with oxygen out of the four fractions. The electrostatic potential (ESP) surfaces for the virgin and aged asphaltenes were determined, which allowed for visualizing the three-dimensional charge distributions of the molecules. The computational method is given in the Appendix A. Figure 13 presents the obtained ESP surfaces, and maximal and minimal ESP points. this purpose since this fraction mostly dominates the engineering properties of asphalt and also presents the highest chemical reactivity with oxygen out of the four fractions. The electrostatic potential (ESP) surfaces for the virgin and aged asphaltenes were determined, which allowed for visualizing the three-dimensional charge distributions of the molecules. The computational method is given in the Appendix A. Figure 13 presents the obtained ESP surfaces, and maximal and minimal ESP points.
As shown in Figure 13a, for the virgin asphaltene molecule, the negative ESP was mostly concentrated at the polyaromatic core due to the π-electrons. The minimum with a value of −21.2 kcal/mol appeared in the vicinity of nitrogen from the pyrrole structure, owing to its lone electron pair. In the aged asphaltene as shown in Figure 13b, the oxygen atoms introduced at the benzyl carbons in the periphery of the polyaromatic core affected the distribution of π-electrons as a result of the strong electronegativity of oxygen and its two lone pairs of electrons. Hence, the local ESP minima were all located around the oxygens, in addition to the pyrrole nitrogen. The global minimum was found around a carbonyl with a value of −35.3 kcal/mol, which was considerably higher than that in the virgin case. The positive ESP was mostly located in the vicinity of the hydrogen atoms. For the virgin asphaltene, the maximal ESP was 33.4 kcal/mol, located along the extension of the N-H bond due to the electronegativity of nitrogen. In the aged case, the maximal ESP appeared at the same location, but was increased to 37.1 kcal/mol. In addition to the ESP surface visualizing the spatial distribution of charges, the dipole moment was also determined for the asphaltene molecules, defined as where μ is the dipole moment vector, and qj and xj are the magnitude and position vector of the j-th charge, respectively. The dipole moment provides a measure of the overall polarity of a molecule. As expected, the aged asphaltene yielded a dipole moment of 5.034 Debye, one order of magnitude higher than 0.419 Debye for the virgin molecule. The above analysis clearly demonstrates that oxidation rendered the charge distribution in the asphaltene molecule much more rugged, extended the positive and negative ESP extremes, and eventually boosted the molecular polarity. Similar observations are expected for the resin and aromatic components. These aging impacts ultimately enhance the electrostatic contribution to molecular interactions within the bulk asphalt and at the asphalt-aggregate interface. As shown in Figure 13a, for the virgin asphaltene molecule, the negative ESP was mostly concentrated at the polyaromatic core due to the π-electrons. The minimum with a value of −21.2 kcal/mol appeared in the vicinity of nitrogen from the pyrrole structure, owing to its lone electron pair. In the aged asphaltene as shown in Figure 13b, the oxygen atoms introduced at the benzyl carbons in the periphery of the polyaromatic core affected the distribution of π-electrons as a result of the strong electronegativity of oxygen and its two lone pairs of electrons. Hence, the local ESP minima were all located around the oxygens, in addition to the pyrrole nitrogen. The global minimum was found around a carbonyl with a value of −35.3 kcal/mol, which was considerably higher than that in the virgin case. The positive ESP was mostly located in the vicinity of the hydrogen atoms. For the virgin asphaltene, the maximal ESP was 33.4 kcal/mol, located along the extension of the N-H bond due to the electronegativity of nitrogen. In the aged case, the maximal ESP appeared at the same location, but was increased to 37.1 kcal/mol.
In addition to the ESP surface visualizing the spatial distribution of charges, the dipole moment was also determined for the asphaltene molecules, defined as where µ is the dipole moment vector, and q j and x j are the magnitude and position vector of the j-th charge, respectively. The dipole moment provides a measure of the overall polarity of a molecule. As expected, the aged asphaltene yielded a dipole moment of 5.034 Debye, one order of magnitude higher than 0.419 Debye for the virgin molecule. The above analysis clearly demonstrates that oxidation rendered the charge distribution in the asphaltene molecule much more rugged, extended the positive and negative ESP extremes, and eventually boosted the molecular polarity. Similar observations are expected for the resin and aromatic components. These aging impacts ultimately enhance the electrostatic contribution to molecular interactions within the bulk asphalt and at the asphalt-aggregate interface.

Conclusions
This study investigated the impacts of oxidative aging on asphalt in terms of the performance-related engineering properties obtained from the molecular dynamics simulation. The following conclusions can be drawn on the basis of the findings:

•
The oxidative aging of asphalt raised the density (from 1.008 to 1.081 g/cm 3 ), cohesive energy density (from 3.02 × 10 8 to 3.49 × 10 8 J/m 3 ), and glass transition temperature (from 315 to 345 K). The overall molecular polarity was elevated upon aging, as reflected in the increased electrostatic contribution to CED. The increase in T g was attributed to the stronger intermolecular attractiveness and reduced molecular mobility, which lowers the resistance of asphalt to thermal cracking. • Aging considerably slowed down the diffusion process (and consequently the selfhealing capacity) of asphalt at a given temperature. The diffusion behaviors of both the virgin and aged asphalt systems at different temperatures were well captured by the Arrhenius relationship, in which prefactor D 0 was a valid index representing the overall effect of aging on asphalt diffusion. Oxidation substantially decreased the prefactor from 1.1 × 10 −5 to 6.8 × 10 −6 cm 2 /s.

•
Aging benefited the adhesion between asphalt and calcite substrate in the dry condition, as it substantially improved the electrostatic interactions at the interface due to the increased molecular polarity. In the wet condition, the use of the proposed index, i.e., the ratio of the residual asphalt-aggregate adhesion to that in the dry condition, indicated that aging yielded consistently higher moisture susceptibility. The simulation results suggest that aging would increase the propensity of asphalt pavement to moisture damage.
The oxygen introduced into the molecules during oxidation significantly affected the charge distributions and ESP extremes. The resulting rise in molecular polarity was highly responsible for the aging-induced changes in the bulk and adhesive properties of asphalt. It can be inferred that, for the purpose of effectively reusing the oxidized asphalts, restoring the aging effects might be approached by mitigating the molecular polarity, which will be explored in future work.
Author Contributions: Conceptualization, methodology, investigation, and writing-original draft preparation, W.C.; methodology, writing-reviewing and editing, E.F. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by the Natural Science Foundation of Hunan Province (CN) (grant no. 2022JJ40589) and Central South University (202045025).

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data that support the findings also form part of an on-going study, and will be available from the corresponding author on request.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The electrostatic potential was calculated on the basis of the Multiwfn code [45,46], and its visualization was accomplished via VMD [47]. The wavefunctions required as input into the code were obtained using Orca 5.0 software [48] with hybrid density functional ωB97M-V [49]. This functional was expected to provide better accuracy than the commonly used B3LYP with dispersion correction (i.e., B3LYP-D3) [50][51][52]. Triple-ζ basis set def2-TZVP was used in calculation [53].