# Model development and prediction of anti-icing longevity of asphalt pavement with salt-storage additive – Journal of Infrastructure Preservation and Resilience

Jan 4, 2022

### Model geometry

Figure 2 depicts the geometries of the overlay in the model. Figure 2a is the cross-sectional image of an asphalt mixture specimen with 5.1 wt% CaCl2-zeolite/p-epoxy #16. The dimension of the cross-section is 150 mm × 32 mm. The geometry of the model was pictured proportionally according to the original cross-sectional image with the aid of the computer-aided design and drafting software AutoCAD, as shown in Fig. 2b. The geometry in Fig. 2b was used to determine the thermal cracking-related coefficient and fatigue cracking-related coefficient of the chloride diffusivity of asphalt mastic. To simplify the geometry, the area in the red box in Fig. 2a and b was adopted for the simulations of water transport and chloride transport, as shown in Fig. 2c. In the Fig. 2b and c, the gray circles are CaCl2-zeolite/p-epoxy particles, the gray irregular polygons are coarse aggregates, and the purple area is the asphalt mastic. Two diameters of the additives were used: 1 mm for the small additives and 2 mm for the big additives. Three asphalt mixture samples were pictured, and the one with the median and thus representative levels of areas of mastic, aggregates, and additives was selected for the modeling study. The software ImageJ was used to measure the areas of mastic, aggregates, and additives in each of the samples, and Table S1 shows the estimated values.

The coarse aggregates in the overlay were assumed to be highly impermeable and thus had little influence on the simulation results of water transport and chloride transport. Therefore, the domain of the coarse aggregates in the geometry was excluded from the simulations.

### Phase transport in porous media (phtr) interface

#### Phase and porous media transport properties of asphalt mastic and additives

In the Brooks and Corey capillary pressure model, the entry capillary pressure (Pec) can be estimated using Eq. 13 [16, 51]:

$${P}_{ec}=frac{2sigma costheta }{{r}_{c}}$$

(13)

where σ is interfacial tension of a fluid, θ is contact angle between the fluid and the capillary tube, and rc is radius of the capillary tube. The σ of water is related to pavement temperatures with a unit of Kelvin (K) [10]. The temperature range for the data in the literature is from 5 °C to 45 °C. The plot of the relationship between the σ of water and pavement temperatures is shown in Figure S1. The contact angles of water to the epoxy layer and the asphalt mastic were measured using water contact angle test. For the pores in the epoxy layer on the surface of the additive, rc is the average value of the pore’s radii measured from scanning electron microscope (SEM) images. For the pores in asphalt mastic, the average rc was calculated with Eq. 14 [41]:

$${r}_{c}=frac{2}{3}*frac{AV%}{100}*{D}_{75}$$

(14)

where AV% is air void percentage or porosity in asphalt mastic, and D75 is the particle-size diameter corresponding to 75% on the cumulative particle-size distribution curve of aggregates (Figure. S2). In the literature, the equation was determined based on the channel theory [38], which is suitable for asphalt pavement with an air void content of higher than 4%. According to Eq. 3, the positive value of the calculated Pec of water in the pores of asphalt mastic was used in the model.

In the Brooks and Corey capillary pressure model, the pore size distribution index of the epoxy layer (λp-epoxy) on the surface of the additive was calculated using Eq. 15:

$${uplambda }_{p-epoxy}=frac{Sigma {(W}_{i}*{A}_{i})}{Sigma {(N}_{i}*{A}_{i})}$$

(15)

where Wi is the weight fraction of pores and Wi = Ni*Ai/∑(Ni*Ai), Ni is the number of pores, and Ai is the surface area of pores. Ni and Ai were obtained from SEM images (e.g., Fig. 3) of the surface of the additive with assisting with AutoCAD. An example of calculating the λp-epoxy for Fig. 3 is shown in Table S2. Three SEM images were analyzed, and the results are shown in Table S3.

The pore size distribution index of asphalt mastic (λp-mastic) was calculated using Eq. 16 [33]:

$${uplambda }_{p-mastic}={(333.3*{k}_{mastic})}^{0.154}$$

(16)

where kmastic is the water permeability of asphalt mastic. Equation 16 was converted from Eq. 2, and the parameters in the equation were determined based on the results of the experiments conducted by the authors of the literature [33].

For simplicity, the density (ρ) of water was approximated to be 1000 kg/m3 for all pavement temperatures. The dynamic viscosity (μ) of water and the ρ and μ of air entailed three equations related to pavement temperatures with a unit of K [27, 42]. In the literature, the temperature range for the μ of water is from 0 °C to 60 °C at atmospheric pressure, and the equations for the ρ and μ of air are valid for a temperature range from -20 °C to 40 °C and a pressure of 1000 bm. The plots of the relationship between the properties and pavement temperatures are shown in Figures S3, S4 and S5. The determination of pavement temperature will be described in the Sect. 3.3.1. The residual saturation (sr) of a water–air system was assumed to be 0 in this study.

The (average) porosity of the epoxy layer on the additive was 2.1%, which can be found in our previous study [53]. The water permeability of the epoxy layer was assumed to be a very small positive value, which was expressed as ‘eps [m2]’ in COMSOL Multiphysics.

The values or equations of the terms mentioned above are shown in Table 1.

#### Initial value and boundary condition

The initial value of the water volume fractions (s0,s1) in the model was assumed to be 0.

The wet-dry cycles of the overlay are mainly affected by precipitations, which are the source of water driving chloride transport. Thus, the hourly data of precipitations including rainfalls and snowfalls from May 1, 2012 to November 27, 2020 was introduced into the model. The weather data at Pullman, Whitman County, Washington was obtained from “AgWeatherNet” of Washington State University (https://weather.wsu.edu/). In the model, the boundary condition on the top surface of the geometry was that either the water volume fraction was 1 for precipitation or the water volume fraction was 0 for no precipitation, which can be implemented using an “Interpolation” function of COMSOL Multiphysics.

### Transport of diluted species (tds) interface

#### Transport properties of asphalt mastic

For a one-dimensional finite diffusion problem, assuming the chloride concentration at surface is constant, the analytical solution of the partial differential Eq. 10 can be expressed as Eq. 17 [9, 19, 22, 23]:

$$frac{Cleft(x,tright)-{C}_{0}}{{C}_{s-{C}_{0}}}=1-mathrm{erf}[frac{x}{2sqrt{Dt}}]$$

(17)

where C(x,t) is the chloride concentration at a location x at a time t, Cs is the chloride concentration at surface, C0 is the initial chloride concentration at location x, erf(z) represents the error function:

$$erf(z) = error function = frac{2}{sqrt{pi }}{int }_{0}^{z}{e}^{{-y}^{2}}dy$$

(18)

Thus, a steady-state chloride diffusion coefficient (Dss) in asphalt mastic with a certain Cs at a given temperature can be estimated with Eqs. 17 and 18 incorporating a mastic diffusion test, which will be described in the Sect. 4.2.

In reality, the D of asphalt pavement is non-steady-state due to the variations of pavement temperature, chloride concentration, the depth, width, and number of thermal cracks and fatigue cracks, which have been considered in the model for more accurate results.

A previous study [37] reported that the D of cement concrete nonlinearly decreases with increasing the square root of C. The plot of the relationship between D and (C + 0.01 mol/m3) of cement concrete with a water-cement ratio of 0.6 is shown in Fig. 4, which was used to represent the relationship between D and C of the studied asphalt overlay, considering their comparable level of chloride diffusivity. A negligibly low background concentration of chloride, 0.01 mol/m3, was introduced to represent chlorides potentially sourced from aggregates and binder, aimed to avoid the non-converged calculation in COMSOL Multiphysics. The water-cement ratio of 0.6 was selected because the D of this cement concrete was the closest to the D of asphalt mixture with 5.1 wt% CaCl2-zeolite/p-epoxy #16 [6, 29, 40, 52]. C-related coefficient to D (Ac) is expressed as:

$${A}_{c} ={(frac{C+0.01}{{C}_{s}})}^{-0.531}$$

(19)

where Cs is the surface chloride concentration used in the mastic diffusion test.

Pavement temperature also has significant influence on D. The effect of pavement temperature on D can be well expressed by the Arrhenius equation [50]:

$$D={D}_{0}bullet mathrm{exp}(-frac{{E}_{a}}{RT})$$

(20)

where D0 is maximal diffusion coefficient, Ea is activation energy for diffusion, R is the universal gas constant, and T is the absolute temperature of asphalt pavement. According to [30], the D and T of cement concrete with a water-cement ratio of 0.6 had the relationship as shown in Fig. 5, which can be used for the relationship between D and T of asphalt pavement. Therefore, T-related coefficient to D (AT) is expressed as:

$${A}_{T} ={e}^{(4008.5times (-frac{1}{T}+frac{1}{{T}_{t}}))}$$

(21)

where Tt is the testing temperature of the mastic diffusion test.

The hourly data of pavement temperature T from May 1, 2012 to November 27, 2020 was calculated with Eq. 22 adapted from [15], and it was implemented using the “Interpolation” function of COMSOL Multiphysics.

$$T=0.686times {T}_{a}+0.000567times {R}_{s}-27.87times d+2.7875$$

(22)

where Ta is air temperature, Rs is daily solar radiation, and d is depth from pavement surface. Arguably the functional pavement might retain more moisture than conventional asphalt pavement and thus feature a slightly different relationship than Eq. (22); but this difference is not considered at this stage of the modeling study.

Moreover, the initiation and propagation of the cracks in cement concrete are known to increase D [8, 12, 45]. In this study, considering the location (Pullman, WA), one can assume that thermal cracking and fatigue cracking were the dominant mechanisms for the development of cracking in the overlay over time. Other risks, such as rutting, raveling, stripping, alligator cracking, were not considered. Similar to cracked cement concrete, the thermal cracks and fatigue cracks in asphalt pavement increase D as well. The thermal cracking mainly results from low temperature environments, while the fatigue cracking is caused by repeated traffic loadings at moderate ambient temperatures.

There is no study reporting how D evolves in asphalt pavement over time as a function of thermal cracking or the relationship between D and the level of thermal cracking. To determine the D-time relationship under the effect of thermal cracking in asphalt pavement, a physical simulation was conducted using COMSOL Multiphysics incorporating a set of thermal cracking-time data adapted from Fig. 13 [14]. Four data points where the thermal cracking percentage of the asphalt pavement was 0% (0 year), 27% (4 years), 38% (9 years), and 59% (11 years) were selected from [14] for this study. According to the definition of the thermal cracking percentage in the reference, the number of the cracks per 150 m length of the survey section was 0, 13.5, 19, and 29.5 corresponding to 0%, 27%, 38%, and 59%, respectively. For a 150 mm length overlay, the number of cracks at 0%, 27%, 38%, and 59% was 0, 0.0135, 0.019, and 0.0295 respectively. In other words, the corresponding number of the 150-mm section per crack was 0, 74, 53, and 34, respectively. To estimate D at different thermal cracking levels, the variation of the average C (Cave) at the bottom of the simulation geometries over time was first determined with Eq. 23:

$${C}_{ave}=frac{{{N}_{n}times C}_{n}+{{N}_{c}times C}_{c}}{{N}_{n}+{N}_{c}}$$

(23)

where Nn is the number of the non-cracking 150-mm section, Cn is the chloride concentration at the bottom of the geometry of a non-cracking 150-mm section, Nc is the number of the 1-crack 150-mm section, and Cc is the chloride concentration at the bottom of the geometry of a 1-crack 150-mm section. Cn and Cc can be obtained from the simulation described below.

Figure 6 shows the geometries used to determine Cn and Cc in COMSOL Multiphysics. The geometry in Fig. 6a is the same as the geometry used in the Sect. 3.1. The geometries in Fig. 6b-d include inverted triangle cracks with different depths and widths. The circles in the geometry were used as traditional aggregates instead of the additive. For simplicity, the model had the following assumptions:

1. 1.

the max crack depth of the overlay was 30 mm at 59%;

2. 2.

the max crack width was half of the max crack depth;

3. 3.

the ratio of two max crack depths was proportional to the ratio of the two thermal cracking percentages.

In the “tds” interface, the aggregates were excluded from computing. The Dss in asphalt mastic had been determined with the mastic diffusion test. It was assumed that the initial value of the chloride concentration in the asphalt mastic was 0. The Cs on the top and crack surfaces of the geometry was 560 mol/m3. The output was the relationship between the chloride concentration (e.g., Cn and Cc) at the bottom of the geometry and time. Eleven points were selected to calculate Cn and Cc over time. The simulated Cave over time at each cracking level was determined with Eq. 23. Meanwhile, the chloride transport follows the Fick’s second law. Thus, D at each cracking level was estimated where the C-time curve calculated according to Eqs. 17 and 18 matched best with the Cave-time curve simulated with COMSOL Multiphysics. Figures S6, S7, S8 and S9 show the simulated and calculated C-time curves and the goodness of fit at different cracking levels. The D-time relationship related to thermal cracking is shown in Fig. 7. Therefore, thermal cracking-related coefficient to D (Atc) is expressed as:

$${A}_{tc}={e}^{(2E-10)times t}$$

(24)

The fatigue cracking-related coefficient to D (Afc) was determined with the same method of Atc. The datapoints of equivalent single axle loads (ESALs) and crack properties (e.g., crack rate and crack width) used for this study were adapted from [21]. Lane #9 was a 90-mm asphalt layer constructed on a 153-mm in-place stabilized soil cement (10%) subbase under a 102-mm layer of crushed stone. The relationship between ESALs and crack rate is shown in Figure S10, which was adapted according to the curve of Lane #9 in Fig. 7.2 in the main text of a reference [21]. The ESALs-crack width relationship as shown in Figure S11 was obtained based on Figures S12 and S13, which were adapted according to the figure of the sample No.25 with a load of 7.25 KN in Appendix D of the reference [21] and the single load equivalency factor (LEF) data of a design guide (Transportation Officials, 1993) respectively. The crack width-crack depth relationship as shown in Figure S14 was adapted from Fig. 5.22 of a reference [21]. The time-ESALs datapoints at State Route 195 near Pullman in Washington State from 2010 to 2019 as shown in Table S4 were obtained from the long-term pavement performance (LTPP) program of Federal Highway Administration, United States Department of Transportation.

It was assumed that the lane width of the 150 mm-section model was 3.66 m and all fatigue cracks were transverse fatigue cracks. According to the equation in Figure S10, Table S4, and Eq. 25, the length of the transverse fatigue crack (L) in a 150 mm-section model was determined. L then was converted to the number of 150-mm sections at different years when one transverse fatigue crack existed and crossed the whole lane of the model, as shown in Table S5. The width of the transverse fatigue crack was determined according to the equation in Figure S11, and the total depth of the transverse fatigue crack was determined according to the equation in Figure S14. It was assumed that the transverse fatigue crack started cracking from the bottom of a 50-mm binder course under the overlay. Thus, the depth of the transverse fatigue crack in the model was calculated using the total depth minus 50 mm. The depth and the width of the transverse fatigue crack are shown in Table S6. In COMSOL Multiphysics, the geometries used to determine Cn and Cc related to the fatigue cracking of asphalt pavement are shown in Fig. 8. The settings in the “tds” interface for fatigue cracking and the calculating procedures were similar with that for thermal cracking.

$$L=pavement surface areatimes fatigue crack rate$$

(25)

Figures S15, S16, S17 and S18 show the simulated and calculated C-time curves and the goodness of fit at different years for fatigue cracking. The D-time relationship related to fatigue cracking is shown in Fig. 9. Therefore, Afc is expressed as:

$${A}_{fc}=1+frac{t}{(8.6E-13)times {(t+e}^{55.5-(1.75E-8)times t})}$$

(26)

Combining Dss, Eq. 19, Eq. 21, Eq. 24, and Eq. 26, the final equation for the D of asphalt mastic (Dmastic) in the tds interface was expressed as:

$${D}_{mastic}={D}_{ss}times {A}_{C}times {A}_{T}times {A}_{tc}times {A}_{fc}$$

(27)

In addition, the threshold of moisture volume fraction (s1) in concrete for chloride diffusion was 20% [35]. As such, Dmastic was expressed as Eq. 27 when s1 was equal or greater than 0.2, while Dmastic was 0 when s1 was less than 0.2, in which s1 was the result of the “phtr” interface.

The parameters related to Dmastic mentioned above are provided in Table 2.

The D in the additives in the presence of water (Dadditive) was estimated with the cooperation of the curve of chloride concentration vs. time in water from an additive immersion test and the curve of chloride concentration vs. time in water simulated with the tds interface of COMSOL Multiphysics. Dadditive was determined when the two curves had the best fit. The additive immersion test is described in the Sect. 4.2.

In the “tds” interface, the geometry used to simulate the C-time curve and the original image of the layer of the compacted asphalt-coated additives in the beaker for the test are shown in Fig. 10a and b, respectively. The total height and the width of the geometry were 132.5 mm and 120 mm respectively, which were the height and the width of the water in the beaker. In the geometry, the gray part is water and the rest represent the compacted asphalt-coated additives with a height of 17.5 mm, which was pictured using AutoCAD according to Fig. 10c. Figure 10c was the processed image according to Fig. 10b using ImageJ.

In the settings of transport properties, Dadditive was a value that needed to be solved to achieve the best match between the experimental vs. simulated chloride concentration curves, and the D in water was 1.2E-9 m2/s derived from the reference [13]. The initial value of the C in water was 0, and the initial value of the C in the additives (Cinitial) was another value to be solved for a match between two chloride concentration curves.

The simulation output was the simulated C-time curve. Simulation results showed the two curves (Figure S19) had the best match when Dadditive was 2E-11 m2/s and Cinitial was 1500 mol/m3. The data points before 71 h in Figure S19b were removed due to the unstable diffusion at the very early stage of the test. Same as the expression of Dmastic, Dadditive was affected by s1 from the result of the “phtr” interface. Dadditive was equal to 2E-11 m2/s when s1 was equal or greater than 0.2, while Dadditive was 0 when s1 was less than 0.2.

#### Initial values and boundary conditions

In reality, CaCl2 exists in a liquid state when water is sufficient, otherwise it is in a solid state. The C in the additives was supposed to be affected by s1 from the water transport model as well, but this effect was difficult to be simulated in COMSOL Multiphysics at this stage. As such, this study assumed that in a very short time the C in the additives reached to the maximum value (7325 mol/m3), which was calculated based on the average volume of the additives. For simplicity, the amount of available water within the sphere of influence of the additives was assumed to be equivalent to the volume of the additives. Subsequently, the C in the additives decreased over time until the system reached a balanced C. To this end, the initial value of the C in the additives in the model was set to be 7325 mol/m3, while the initial value of the C in the asphalt mastic was assumed to be 0 mol/m3.

The C in asphalt mastic and on the top surface of the overlay may gradually increase over time until it reaches to the maximum C. In the service environment, part of the salt on overlay surface is washed away by rainfalls and snowfalls. There is no study reporting how much of the salt is removed by rainfalls or snowfalls. Based on field observations, this study assumes that approximate 5% and 50% of the C on the top surface of the overlay remained after the hour of rainfall and snowfall, respectively. The boundary condition related to C on the top surface of the geometry was implemented by the functions of “Flux” and “Interpolation” in COMSOL Multiphysics. The snowfall season was defined as from November 13 to February 25 every year in Pullman, Washington.

Figure 11 shows the flowchart of the water transport model, the chloride transport model, and how the output of the water transport model affects the chloride transport model.

### Mesh and study

The asphalt mastic, additives, and top surface boundary were meshed according to their sizes. The free triangular mesh was selected as the mesh of the model. This study selected the normal size, fine size, and finer size of mesh for the asphalt mastic, additives, and top surface boundary, respectively.

Both the water transport time-dependent model and the chloride transport time-dependent model were studied for 44,280 h (approximate 5.05 years) with an interval of 1 h. This is because in practice such an overlay would be replaced after 5 years in service and some of the model assumptions would no longer hold true if the asphalt overlay serves beyond 5 years in the field where moisture damage, UV aging, etc. start to notably compromise the integrity of the overlay.