Quantifying oxygen diffusion in bitumen films using molecular dynamics simulations

.


Introduction
Bitumen, a residue from the distillation of crude oil, is widely used as a binder in asphalt mixtures for paving road surfaces.The bitumen's durability plays a vital role in the mechanical property of asphalt mixtures thus can significantly influence asphalt pavement's service life.The oxidative ageing of the bitumen has been one of the major concerns, which can lead to severe hardening and brittleness of asphalt mixtures and eventually cause the distresses such as raveling and block cracking of asphalt pavements [1][2][3][4].
Researchers have investigated the effects of oxidative ageing on the rheological and viscoelastic properties of bituminous materials.Camargo et al. [5] presented a study on the effects of oxidative aging on the rheology and chemical composition of bitumen using four laboratory aging procedures, i.e., rolling thin-film oven (RTFO), pressure aging vessel (PAV), nitrogen atmosphere oven aging test (NAAT) and ambient atmosphere oven aging test (OAAT).It was found that the oxidative aging was the main cause for the bitumen's hardening.Based on the Fourier Transform infrared spectroscopy (FTIR) and Dynamic Shear Rheometer (DSR) tests, Liu et al. [6] developed the chemical-rheological property relationship of bitumen subjected to the oxidative aging.The carbonyl index was well linked to the complex shear modulus and phase angle of the bitumen.Alavi et al. [7] proposed an approach to quantify the impact of bitumen oxidative aging on the viscoelastic properties of asphalt mixtures.The complex modulus of the asphalt mixture and the carbonyl area of the recovered bitumen were measured after the mixtures were subjected to the laboratory aging with the varying durations.A good correlation was found between the carbonyl in the bitumen and the parameters of continuous relaxation spectrum for the asphalt mixture.
The mechanisms of oxidative ageing have been explored in a number of existing studies [3,8,9].The oxidative ageing is an irreversible diffusion-driven phenomenon mainly controlled by chemical reactions between bitumen components and atmospheric oxygen when the surface of asphalt pavement is exposed to the atmosphere.The oxidation process alters the chemical features of the bitumen, e.g., the formation of carbonyl and sulfoxide functional groups [3], which has a significant influence on the physical and mechanical properties of bituminous materials including rheology [10], adhesion [11,12] as well as fatigue [13][14][15].During the whole ageing process, oxygen molecules need to first diffuse into the bitumen film to a certain depth and then react with the bitumen molecules at that position, thus the oxidative reaction is strongly dependent on the physical diffusion of the oxygen into the bitumen film [1,[16][17][18].Furthermore, researchers have been trying to employ the anti-ageing compounds (AACs) to inhibit the bitumen's oxidative ageing [19][20][21].To select the potential AACs, a screening method was proposed by Omairey et al. [20] to evaluate the anti-ageing effectiveness of different types of AACs.However, the anti-ageing mechanism of AACs is still not known, especially at the molecular scale.
The oxygen diffusivity value in the bitumen is a crucial parameter in the transport models that are used for estimating bitumen oxidation in pavements [22,23].The literature is rather limited on the diffusion coefficient of oxygen in the bitumen.The difficulty in measuring the diffusion coefficient is separating physical diffusion behaviour from the chemical reaction during the complicated diffusion-reaction process of bitumen oxidative ageing.Perisamy [24] directly measured the oxygen diffusion coefficients in several bitumens using mass changes, but the oxidative reaction with the bitumen was not considered in the measurement.To avoid the effect of the reaction, Dickinson [25] and Jing et al. [26] used nitrogen to replace oxygen in the evaluation of the oxygen diffusion coefficient in the bitumen.Glover et al. [23,27] indirectly determined the oxygen diffusivity based on laboratory oxidation experiments in bitumen films of known reaction kinetics.Therefore, it is still extremely challenging to characterise the physical diffusion of oxygen through the experimental methods due to the chemical interactions between oxygen and bitumen.While the oxygen molecules diffuse into the bitumen film, they can be consumed by the oxidative reaction with the bitumen molecules, which masks the real diffusion process and thus significantly affects the accuracy of the determination of the oxygen diffusion coefficient.
With the recent advances in the high performance computation, molecular dynamics (MD) simulation has become a powerful tool to investigate the thermodynamic behaviours of materials at the molecular scale [28][29][30], which has a strong potential to provide a viable solution for the accurate determination of the oxygen diffusion coefficient in the bitumen.Tsige et al. [31] used the MD simulation to study the diffusion of a solvent into a polymer melt based on Fickian diffusion.It was found that the diffusivity was essentially constant for the low and intermediate solvent concentration.The MD simulation was also applied by Neyertz et al. [32] and Lin et al. [33] to calculate the diffusion coefficients of small molecules into the polymer films.To understand the decomposition of covalent adaptable network, Sun et al. [34] employed the MD simulation to investigate the diffusion of ethylene glycol (EG) solvent into the polymer film composed of DGEBA (Diglycidyl ether of bisphenol A) and fatty acid.Recently, research efforts have primarily used the MD simulation to investigate the diffusion between virgin bitumen [35,36] or rejuvenator [37,38] and aged bitumen as well as the moisture diffusion in the bituminous materials [39,40].These studies indicate that the MD simulation can be used to fundamentally interpret the diffusion behaviours of molecules in the materials.
The objective of this study is to use the MD simulation to investigate the physical diffusion of the oxygen into the bitumen film and antiageing compounds' effects on the oxygen diffusion.The molecular model of the bitumen film was firstly constructed with the saturate, aromatic, resin and asphaltene (SARA) four fractions.A bitumen-air bilayer model was then developed by combining the bitumen film with air models.The MD simulations were performed for 5 ns (nanosecond) at different temperatures (i.e., 25, 50 and 100 • C) to simulate the physical diffusion process of the oxygen into the bitumen film.Based on the simulation results, Fick's second law was used to calculate the diffusion coefficient of the oxygen in the bitumen film.Furthermore, the effects of the temperature and anti-ageing compounds were analysed based on the proposed model.Fig. 1 shows an overview of the research steps taken in this study.

Molecular models of bitumen
Bitumen is a viscoelastic material with a very complex chemical composition.Based on the physical and solubility properties, the bitumen is generally divided into saturates (S), aromatics (A), resins (R) and asphaltenes (A) four fractions.A 12-component bitumen molecular model has been proposed by Li and Greenfield [41] to represent Strategic Highway Research Program (SHRP) core bitumen AAA-1 with a penetration of 150/200.Fig. 2 illustrates the molecular structures of the twelve components in this model.The asphaltenes (A) contain asphaltene-phenol, asphaltene-pyrrole and asphaltene-thiophene. Five different types of molecules (i.e., quinolinohopane, thioisorenieratane, benzobisbenzothiophene, pyridinohopane and trimethylbenzeneoxane) are selected to represent the resins (R).The aromatics (A) are composed of perhydrophenanthrene-naphthalene (PHPN) and dioctylcyclohexane-naphthalene (DOCHN).Squalene and hopane are used for the saturates (S).Table 1 shows the chemical composition of the bitumen molecular model.
The twelve different types of molecules shown in Fig. 2 was used in this study to construct the bitumen model.Based on their numbers in Table 1, these molecules were put in a 3D periodical simulation box with a low initial density of 0.1 g/cm 3 to form a bulk bitumen model.A geometry optimisation was run for 5000 iterations to minimize the system energy and obtain a better initial configuration.The system was then subjected to a molecular dynamic (MD) equilibration under an NVT ensemble (constant number of particles N, constant system volume V and constant temperature T) at 298 K for 100 ps (picosecond).A MD simulation with an NPT ensemble (constant pressure P) was followed to relax and refine the configuration for 200 ps at 1 atm.The simulation time is sufficient as both the density and energy of the system have approached to a stable state, as shown in Fig. 3

Molecular models of bitumen-air bi-layer system
In addition to the bitumen molecular model, an air molecular model is needed for the diffusion simulations of the oxygen in the atmosphere into the bitumen film.In reality, the air is composed of 78 % nitrogen and 21 % oxygen with small amounts of other gases (such as carbon dioxide, neon, and hydrogen).In this study, a 2-component air molecular model was used to represent the air, which was made up of the nitrogen and oxygen molecules.
To construct a bitumen-air bi-layer model, two confined layers were first prepared for the air and the bitumen, respectively.The confined layer had a flat surface in the z direction, which was different from the bulk model.A confined bitumen layer with a certain density that was the same as that of the equilibrated bulk model obtained in Section 2.1 was established by randomly assembly molecules into a simulation box.A confined air layer containing 3000 nitrogen molecules and 800 oxygen molecules was prepared as well.The air density in the simulation was higher than the real air density to ensure that the system maintains a sufficient oxygen supply for diffusion in the bitumen film at the small scale.It should be noted that the diffusion coefficient is usually independent of the concentration [37,38,42].The higher oxygen levels than the actual values have been used in the molecular simulations of bitumen oxidative ageing [43,44].The initial confined layers were first subjected to a geometry optimization to minimize system energy.A MD simulation with the NVT ensemble was then conducted at 298 K for 200 ps to equilibrate the systems.The equilibrated models of the air and bitumen layers are shown in Fig. 4(a) and (b), respectively.
Once the confined layers were obtained, a bitumen-air bi-layered system model can be built by attaching the air layer to the bitumen layer.The bitumen-air bi-layer model with a dimension of 39 Å (angstrom) × 39 Å × 163 Å is illustrated in Fig. 4(c).

Diffusion simulations
Molecular dynamics (MD) simulations were performed on the bitumen-air bi-layer system model developed in Section 2.2 to simulate the diffusion of the oxygen into the bitumen film.The bi-layer model shown in Fig. 4(c) had a periodic boundary condition in the x and y directions, while the non-periodic condition was applied in the z direction (diffusion direction).Two reflective walls were added at the top and bottom edges of the system in the z direction to prevent atoms to move outside the edges.It is noted that the interactions between the atoms in the system and the walls are very weak and thus have negligible effect on the following MD simulations.A geometry optimization was first conducted on the bi-layer model.After that, the model was equilibrated by a MD simulation in the NVT ensemble.Finally, an additional MD simulation with the NVT ensemble was performed for 5 ns to simulate the diffusion process.The diffusion simulation time (5 ns) was determined by six trail simulations.The diffusion coefficients of oxygen were calculated for different simulation times (1, 2, 3, 4, 5 and 6 ns).It was found that the diffusion coefficients for 4, 5 and 6 ns are nearly the same, which means that the diffusion coefficient trends to be constant after 4 ns.In order to save simulation time, 5 ns was used for diffusion simulation in this study.
In this study, all the molecular models were prepared using Materials Studio [45] and then transferred to Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) for the subsequent MD simulations [46].The atomic trajectories of the diffusion process were viewed through the visualization tool OVITO [47].All the simulations  were performed with a time step of 1 fs (femtosecond).Atom-based summation with a cut-off distance of 15.5 Å was applied for van der Waals interactions between molecules, and Ewald summation method was employed for the electrostatic interactions.The temperature and pressure of the simulation system were controlled by Nose-Hoover-Langevin (NHL) thermostat and Andersen barostat, respectively.

Interatomic potentials (force field)
In the simulation systems, the interactions between the atoms are described by a Polymer Consistent Force Field (PCFF) that has been extensively parameterized and validated for application to polymers and organic materials [48][49][50][51][52].The total potential energy (E total ) in the PCFF is composed of both valence and non-bond terms, as shown in Eq. (1).
The valence terms contain the interactions caused by bond (E b ), angle (E a ), torsion angle (E t ) and out-of-plane angle (E o ), as well as coupling between them (E bb ′ , E ba , E bt , E aa ′ , E at and E aa ′ t ).The non-bond terms consist of the van der Waals forces (E vdw ) represented by the Lennard-Jones (LJ) 9-6 potential and the Coulombic forces (E elec ) between partially charged atoms, as expressed in Eqs. ( 2) and (3).
where ε ij is the coefficient of well-depth energy that represents the strength of the integrations, r ij is the distance between atoms i and j, r 0 ij is the distance at which the LJ potential is zero, and q i and q j are the partial charges on atoms i and j, respectively.
It should be noted that the PCFF only defines the physical interactions between the atoms without the chemical reactions.Therefore, using the non-reactive force field, the diffusion simulation in this study is a physical process, which means that the oxidative reactions will not occur in the bitumen-air bi-layer system during the whole simulation.This enables one to separate physical diffusion behaviour from the chemical reaction, which is extremely difficult for the experimental measurement of the oxygen diffusion coefficient in the bitumen.

Observation of diffusion process
To directly observe the diffusion process of oxygen molecules into the bitumen film, the molecular trajectories of the bitumen-air bi-layer system were recorded during the MD simulations.Fig. 5 shows the snapshots of the bi-layer system at the simulation time 0 ns, 1.5 ns and 5 ns at 50 • C. In the simulation model, the upper area is the air composed of the nitrogen (in green) and oxygen (in red) molecules and the lower area is the bitumen film with 12 components.It can be seen form Fig. 5 (a) that the contact interface (indicated by the dash line) between the air and the bitumen film is relatively smooth at the beginning of simulation (t = 0 ns).After the simulation of 1.5 ns, the air molecules have moved across the original interface and penetrated the bitumen film, as shown in Fig. 5(b).To better display the oxygen diffusion, the bitumen and nitrogen molecules were set to be invisible in Fig. 5(c) where only oxygen molecules were shown.It can be clearly observed that the oxygen molecules have moved into the half of the bitumen area.Fig. 5(d) and (e) present the microstructures of the bitumen-air system after 5 ns.It is found that the oxygen molecules further diffused and have penetrated the whole bitumen film.The air molecules immersed much deeper into the bitumen film than the bitumen molecules into the air layer, which shows asymmetric diffusion characteristic.This is reasonable because the air has smaller molecular size and thus has higher mobility than the bitumen.
To quantify the diffusion, the time-dependent mass density profiles for oxygen molecules in the bitumen-air bi-layer system were monitored  during the MD simulation based on the one-dimensional (1D) binning analysis along the diffusion (z) direction with a thickness of 3 Å.Fig. 6 shows the evolution of mass density profiles of the oxygen molecules.The density profiles at different simulation time (1 ns, 2 ns, 3 ns, 4 ns and 5 ns) have similar shapes.The mass density of oxygen gradually increased at the lower z positions (bitumen-rich side) with the simulation time, which demonstrates that the oxygen molecules are diffusing into the bitumen film.The penetration depth can be calculated as the difference between the z positions at different simulation time.

Calculation of diffusion coefficient
Diffusion is the molecular movement from a region of high concentration to a region of low concentration.It is one-dimensional, driven by a gradient in the concentration.Fick's second law [53] shown in Eq. ( 4) is often used to describe the diffusion.It is assumed in this equation that the diffusion does not change the volume of the medium.The nature of the diffusion is called Fickian in this case.
where c is concentration; t is time; z is position; and D(c) is diffusion coefficient.
In general, the diffusion coefficient D(c) is assumed to be independent of the concentration c and considered as a constant D 0 , thus Eq. ( 4) can be expressed as a function of time and position, as shown in Eq. ( 5).In this study, the diffusion coefficient is treated as a constant and will be calculated by fitting expression in Eq. ( 5) to the MD simulation results based on the bitumen-air bi-layer model.The functional form Eq. ( 5) has been widely applied for polymers [33,42,54] and bituminous materials [37,38] to calculate their diffusion coefficients.
where c 0 is the equilibrium oxygen concentration; erf is the error function; and D 0 is the constant diffusion coefficient.
Here Eq. ( 5) was used to calculate the diffusion coefficient of oxygen into the bitumen film.It is clear from Eq. ( 5) that the diffusion coefficient can be extracted from the mass density profiles scaled with zt − 0.5 .Fig. 7 illustrates the mass density profiles scaled with zt − 0.5 for oxygen molecules in the bitumen-air bi-layer system.It can be seen that the density profiles at different simulation time form a single master curve, which indicates that the oxygen diffusion into the bitumen film follows the Fickian diffusion behaviours [31].
The diffusion coefficient was obtained by fitting the scaled density data from MD simulations using Eq. ( 5), as shown in Fig. 7.The calculated diffusion coefficient at 50 • C is 2.01 × 10 − 10 m 2 /s for the oxygen in the bitumen-air bi-layer system.The result agrees well with the findings from the experimental study [16] in which the oxygen diffusion coefficient into the PEN 180/200 bitumen was 1.83 × 10 − 10 m 2 /s.It is noted that the penetration of PEN 180/200 bitumen is really close to that (150/200) of the bitumen AAA-1.Han et al. [27] also reported that the oxygen diffusivities ranged from 10 − 10 to 10 − 11 m 2 /s for the neat bitumens with different PG Grades (five base binders and three polymermodified binders).This demonstrates that the MD simulating approaches is capable of reliably determining the oxygen diffusion coefficient by avoid the time-consuming and costly experimental tests.It can also be employed, from the molecular perspective, to understand the physical diffusion mechanism of oxygen in the bitumen film.
The diffusion behaviour of nitrogen into bitumen was also analysed following the same simulation procedure.The calculated diffusion coefficient of nitrogen is 1.63 × 10 − 10 m 2 /s at 50 • C. It is seen that nitrogen has smaller diffusion coefficient than oxygen (2.01 × 10 − 10 m 2 /s), which indicates that the nitrogen diffuses slower in the bitumen.

Effect of temperature on diffusion
The MD diffusion simulations were also performed at different temperatures to investigate the effect of temperature on the diffusion of oxygen molecules into the bitumen film.Fig. 8 shows the positions of oxygen molecules in the bitumen-air bi-layer model after the diffusion simulation of 5 ns at the temperatures of 25, 50 and 100 • C. To better display the oxygen diffusion, it is noted that only oxygen molecules were shown in the bitumen-air bi-layer model in Fig. 8 where the bitumen and nitrogen molecules were set to be invisible.Compared to the initial position of oxygen molecules in the model shown in Fig. 8(a), more the oxygen molecules have moved to the lower area (the bitumen area) of   the simulation box at all three temperatures (see Fig. 8(b)-(d)) after the simulation of 5 ns.It is clearly observed that the amount of oxygen molecules in the bitumen area increased with the increasing temperature, which was quantified by the mass density profiles of oxygen molecules at simulation time of 5 ns shown in Fig. 9. Based on the MD simulation results, the diffusion coefficients of oxygen molecules in the bitumen-air bi-layer system were calculated at different temperatures using Eq. ( 5) and shown in Fig. 10.The diffusion coefficients of oxygen molecules at the temperatures of 25, 50 and 100 • C were 7.68 × 10 − 11 m 2 /s, 2.01 × 10 − 10 m 2 /s and 5.81 × 10 − 10 m 2 /s, respectively.It can be seen that the diffusion coefficients increased as the temperature increased, which indicates that the diffusion of oxygen molecules into the bitumen film is positively temperature dependent.This finding is consistent with laboratory measurements, in which Glover et al. [23] have reported that the oxygen diffusion coefficient is a function of the temperature.The higher the temperature, the faster the oxygen diffusion into the bitumen film.This is because that the oxygen molecules have the more kinetic energy and thus move more quickly at a higher temperature.

Impact of anti-ageing compounds (AACs) on oxygen diffusion
Based on the proposed MD method, this section will investigate the impact of AACs on the physical diffusion of oxygen molecules into the bitumen film to explore the anti-ageing mechanism of AACs.

Selection of AACs
In this study, two types of AACs, Irganox acid and Dilauryl thiodipropionate (DLTDP):furfural, were selected based on the authors' previous work [10].Irganox acid synthesised in the laboratory [55] is a hindered phenol, which is typically used as a precursor for polymer antioxidants.DLTDP:furfural is the commercial compounds.Furfural is an organic aldehyde and DLTDP is generally applied as a secondary antioxidant in polymers.Table 2 shows the physicochemical characteristics of the three additives.
Gao et al. [10] have evaluated the anti-ageing effectiveness of Irganox acid and DLTDP:furfural using Fourier Transformation Infrared (FTIR) spectroscopy.In their studies, the concentration of 10% (by weight of the base bitumen) was used for Irganox acid to fabricate the AAC-modified bitumen.The 1.5% DLTDP and 2% furfural were employed for the DLTDP:furfural formulation.It was found that both Irganox acid and DLTDP:furfural exhibited a high anti-ageing performance.Based on the concentration of AACs used for fabricating the AAC-modified bitumen in the previous experimental work [10], the number of the AAC molecules can be determined for the MD simulations, as shown in Table 3.These AAC molecules together with the bitumen molecules shown in Fig. 2 were put in a simulation box to build the AAC-modified bitumen model.Following the same procedure as the bitumen model and the bitumen-air bi-layer model described in Sections 2 and 3, the AAC-modified bitumen model and the AAC-modified bitumen-air bi-layer model were constructed and simulated to investigate the physical diffusion of oxygen into the AAC-modified bitumen film.Fig. 12 shows the initial Irganox acid-modified bitumen-air bi-layer model.

Diffusion of oxygen into AAC-modified bitumen film
After the MD simulations, the diffusion coefficients of oxygen molecules into the AAC-modified bitumen-air bi-layer system were calculated at different temperatures using Eq. ( 5), as shown in Fig. 13.It is found from Fig. 13 that the oxygen diffusion coefficients for both the Irganox acid and DLTDP:furfural modified bitumen samples significantly increased with an increase in the temperature, which is consistent with that for the reference bitumen.The oxygen diffused faster at a higher temperature due to more kinetic energy.This further proves that the diffusion of oxygen molecules into the bitumen film is strongly temperature dependent.
The diffusion coefficients of oxygen molecules into the Irganox acidmodified bitumen film at the temperatures of 25, 50 and 100 • C were 7.45 × 10 − 11 m 2 /s, 1.43 × 10 − 10 m 2 /s and 3.64 × 10 − 10 m 2 /s, respectively.It can be clearly seen from Fig. 13 that the Irganox acidmodified bitumen has smaller diffusion coefficients than the reference   bitumen, especially at the higher temperatures (e.g., 50 and 100 • C).This demonstrates that the oxygen diffused slower in the Irganox acidmodified bitumen.Irganox acid reduced the oxygen diffusion, leading to the lower oxygen concentration in the bitumen film.It seriously limited the oxidative reaction by preventing oxygen to reach the freeradicals of the bitumen film and ultimately hindered the ageing of the bitumen.Thus, it is concluded that reducing the oxygen physical diffusion into the bitumen film is one of the anti-ageing mechanisms of Irganox acid.
For DLTDP:furfural-modified bitumen, the diffusion coefficients of oxygen molecules were 1.13 × 10 − 10 m 2 /s, 2.38 × 10 − 10 m 2 /s and 6.67 × 10 − 10 m 2 /s at the temperatures of 25, 50 and 100 • C, respectively.Compared to the reference bitumen, all the oxygen diffusion coefficients are a little higher for the DLTDP:furfural-modified bitumen, which indicates that DLTDP:furfural slightly accelerated the oxygen diffusion in the bitumen film.Thus, the physical diffusion is not a contributor to the anti-ageing behaviours of DLTDP:furfural.In fact, controlling the oxidative reaction between the oxgen and the bitumen film is an important anti-ageing mechanism of DLTDP:furfural.It has been reported by Omairey et al. [56] that in the antioxidant treatment furfural is a primary antioxidant used as a free-radical scavenger, and DLTDP is a secondary antioxidant working as a peroxide inhibitor.The treatment containing a combination of furfural and DLTDP involves a two-stage reaction, i.e., condensation reaction of furfural and bitumen and peroxide decomposition reaction of DLTDP and bitumen [19].

Inhibition of bitumen aging by reducing oxygen physical diffusion
As discussed above, Irganox acid and DLTDP:furfural show two different anti-aging mechanisms including 1) reducing the oxygen physical diffusion and 2) controlling the chemical oxidative reaction.Based on the second mechanism, a widely used method to inhibit the bitumen aging is adding antioxidants into the bitumen where the antioxidants react with in-process oxidation product.However, the limitations of this method are: 1) the antioxidants are not applicable for all the bitumens; and 2) the protection function of this kind of antioxidants would be lost as they were gradually depleted.
From the first mechanism, a new method is proposed in this study, which is reducing the oxygen diffusivity by constructing a network in the bitumen to retard the oxygen penetration and increase the transport path, as that has functioned by Irganox acid.This proposed method overcomes the limitations of the existing method, which can be used for all types of bitumens without the antioxidant consumption.Thus, the MD computational framework established in this study can contribute to instruct developing new antioxidant in addition to evaluating the oxygen physical diffusion of the existing bitumen.

Conclusions
In this study, a physical diffusion of oxygen into bitumen film was investigated using molecular dynamics (MD) simulations.A bitumen-air bi-layer molecular model was firstly constructed.The MD simulations with the Polymer Consistent Force Field (PCFF) were then performed to study the oxygen diffusion process in the unmodified bitumen and modified bitumen with two different anti-ageing compounds (AACs) at different temperatures.The diffusion coefficient of the oxygen in the bitumen film was calculated based on the diffusion theory and simulation data.The developed method was also used to analyse the effects of AACs on the oxygen diffusion in the bitumen.The following conclusions can be drawn from this research: (1) The proposed computational approach based on MD simulation is capable of determining the oxygen diffusion coefficient in the bitumen by separating physical diffusion from the chemical oxidative reaction, which is extremely challenging for the experimental measurement but of technological importance.(2) The oxygen diffusion coefficients ranged from 6.67 × 10 − 10 to 7.45 × 10 − 11 m 2 /s for the unmodified and AAC-modified bitumens at the testing temperatures of 25, 50 and 100 • C. (3) The higher the temperature, the faster the oxygen diffusion for the unmodified and AAC-modified bitumens.This proves that the oxygen diffusivity is strongly temperature dependent.(4) The bitumen modified by Irganox acid had a smaller coefficient of the oxygen diffusion, which means that the Irganox acid effectively hindered the oxygen diffusion into the bitumen film.
Reducing the oxygen physical diffusion into the bitumen film is one of the anti-ageing mechanisms of Irganox acid.(5) The diffusion coefficient becomes higher when the bitumen is modified by the DLTDP:furfural.Thus the oxygen physical diffusion is not a contributor to the anti-ageing behaviours of DLTDP:furfural.(6) Reducing the oxygen diffusivity by constructing a network in the bitumen to retard oxygen penetration and increase the transport path is an efficient way to slow down the bitumen aging.It can overcome the limitations of the existing method, which can be used for all types of bitumens without the antioxidant consumption.The MD-based computational approach established in this study can contribute to instruct developing new antioxidant.
Future study will be focused on performing more oxygen diffusion simulations for different bitumens.With the development of advanced experimental techniques, their oxygen diffusion coefficients also will be measured in the laboratory to further validate the proposed MD simulation approach.

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.

Fig. 1 .
Fig. 1.Flowchart of research steps in this study.

Fig. 2 .
Fig. 2. Molecular structures of 12 components in bitumen model.Hydrogen atoms are white, carbon atoms are grey, oxygen atoms are red, nitrogen atoms are blue, and yellow is sulphur atoms.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Y
.Gao et al.

Fig. 6 .
Fig. 6.Mass density profile of oxygen molecules during diffusion at different diffusion times and 50 • C.

Fig. 11
Fig. 11 presents the molecular models of three additives used in the AACs.The chemical formula of Irganox acid is C 17 H 26 O 3 , including one benzene ring.Furfural (C 5 H 4 O 2 ) has a small molecular size and weight.The molecular structure of DLTDP (C 30 H 58 O 4 S) is symmetric with two long alkane (C 12 H 25 ) function groups.Based on the concentration of AACs used for fabricating the AAC-modified bitumen in the previous experimental work[10], the number of the AAC molecules can be determined for the MD simulations, as shown in Table3.These AAC molecules together with the bitumen molecules shown in Fig.2were put in a simulation box to build the AAC-modified bitumen model.Following the same procedure as the bitumen model and the bitumen-air bi-layer model described in Sections 2 and 3, the AAC-modified bitumen model and the AAC-modified bitumen-air bi-layer model were constructed and simulated to investigate the physical diffusion of oxygen into the AAC-modified bitumen film.Fig.12shows the initial Irganox acid-modified bitumen-air bi-layer model.

Fig. 11 .
Fig. 11.Molecular structures of three additives used in the AACs.

Fig. 12 .Fig. 13 .
Fig. 12. Irganox acid-modified bitumen-air bi-layer model for MD diffusion simulation.Molecules of Irganox acid are shown in yellow.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1
Chemical composition of bitumen molecular model used for MD simulations.
Y.Gao et al.

Table 2
Physicochemical characteristics of three additives used in AACs.

Table 3
Chemical composition of AACs used for construction of AAC-modified bitumen model.