Model development and prediction of anti-icing longevity of asphalt pavement with salt-storage additive

This study established a systematic simulation framework to predict the anti-icing longevity of a thin overlay of asphalt pavement with salt-storage additive (APSSA). The water and chloride transport in the overlay when subjected to varying precipitation, temperature, thermal cracking, and fatigue cracking over time were modeled using a Finite Element Method based software. The simulation included two parts: water transport followed by chloride transport. Water transport that obeys the law of conservation of mass was modeled using the phase transport in porous media (phtr) interface of COMSOL, while chloride transport based on Fick’s second law was modeled with the transport of diluted species (tds) interface. The simulation results show that the anti-icing function of a 16-mm thick overlay was fully effective in 2 years and 5 years for the minimum pavement temperature above -3.4 °C and -2.4 °C, respectively. These two pavement temperatures are equivalent to 97.4-percentile and 96.3-percentile of historical hourly pavement temperature near Pullman, Washington.


Introduction
A few studies have reported that the anti-icing life of asphalt pavement with salt-storage additive (APSSA) can be calculated with predictive equations incorporating laboratory experiments [11,24,25,[55][56][57][58]. However, these estimations were overly simplified, rough, and inaccurate because of lack of a mechanistic model. Another significant deficiency of those predictions is neglect of a comprehensive consideration of the factors influencing the anti-icing life of APSSA, such as environmental conditions (e.g., precipitation and temperature) and pavement conditions (e.g., thermal cracking and fatigue cracking). Although many studies have worked on modeling ions (e.g., chloride) transport from the outside to the inside of cement concrete [2,9,18,22,31,36,49], no one has simulated outward ions (e.g., chloride) transport in asphalt pavement or cement concrete. No model has been developed to predict the anti-icing life of APSSA. Therefore, it is necessary to establish a systematic model for a more mechanistic and reliable prediction of the anti-icing life of APSSA.
The COMSOL Multiphysics is a cross-platform finite element method (FEM) based multiphysics simulation software. This software supports the simulation of designs involving single or coupled physics such as electromagnetics, structural mechanics, acoustics, fluid flow, heat transfer, chemical engineering, etc. [7,17,20,26,32,36,47]. Among them, the phase transport in porous media (phtr) interface in the subsurface flow module and the transport of diluted species (tds) interface in the chemical reaction engineering module have been widely applied for simulating multi-phase flow (e.g., water-air flow) in porous media [28,46,48,51] and simulating ionic ingress into cement concrete due to concentration gradient [2,9,22,36,49], respectively. For APSSA, ionic transport from the inside to the surface must be facilitated by water, and the process include two activities: water transport in asphalt mastic and salt-storage additive (SSA) followed by ion (e.g., chloride) transport due to concentration gradient. The two activities can be well modeled and linked by the cooperation of the two interfaces of the FEM software.
This work pioneered a numerical model with the COM-SOL Multiphysics to predict the anti-icing life of a thin overlay of APSSA when subjected to varying precipitation, temperature, thermal cracking, and fatigue cracking over time. In this study, the APSSA sample used for the experiments to obtain the parameters of the model was asphalt mixture with 5.1 wt% laboratory-prepared functional additive (CaCl 2 -zeolite/p-epoxy) #16. The additive was fabricated through a surface treatment approach, in which zeolite particles containing calcium chloride were coated by a microporous epoxy layer. The #16 denotes the gradation of the additive, which passed the sieve #16 and stayed on the sieve #30 (1.2 mm -2.4 mm). More details of the preparation of CaCl 2 -zeolite/p-epoxy #16 and the APSSA sample can be found in the laboratory study [53]. The simulation was separated into two steps, which are water transport into the overlay first then followed by chloride transport from the inside to the surface. The Graphical abstract Page 3 of 28 Zhang et al. J Infrastruct Preserv Resil (2022) 3:1 water transport process follows the law of conservation of mass, while the chloride transport from high-concentration areas to low-concentration areas obeys Fick's second law.

Model description
To activate the anti-icing function of the functional additives in the overlay, there are two critical processes: 1) wetting of solid calcium chloride (CaCl 2 ) in the additives and 2) outward transport of the chloride ions. The wetting process is the one by which water intrudes into the asphalt mastic from the surface and finally enters the functional additives first through various pathways (e.g., pores and defects) in the overlay and then through the micropores on the epoxy coating layer. Subsequently, CaCl 2 turns from its combined state (i.e., solid CaCl 2 inside the additives) to ionic state (Ca cations and Cl anions dissolved in water) in the additives. Thereafter, the chloride ions in water migrate to the surface of the overlay through the same pathways of water ingress, mainly by diffusion (due to the chloride concentration gradient inside the overlay). The schematic illustration of the two processes is shown in Fig. 1.
In reality, the overlay counters wet-dry cycles due to weather variations over time. The wet-dry cycles lead to water-air transport in the overlay, which is a problem of a two-phase flow in porous media that obeys the law of conservation of mass. On the other hand, Fick's second law is used to describe the transporting process of the chloride ions in the overlay, which is a non-steady state diffusion due to the chloride concentration gradient between the inside and the surface of the overlay.
Wet-dry cycles of the overlay Sufficient water in the overlay is an essential precondition that makes the chloride ions being transported outward from inside of the additives to the surface of the overlay. The sources of water are rainfall in summer and snowfall in winter. To solve water-air phase transport in the overlay, the macroscopic mass conservation equation was used [3,43,44,46]: where ε p is the porosity of porous media, ρ i is the density of phase i, s i is the volume fraction of phase i, and u i is the volumetric flux (velocity vector) of phase i. u i is determined with the extended Darcy's law as shown in Eq. 2 [5,44]:  where k ri is the relative permeability of phase i, u i is the dynamic viscosity of phase i, k is the permeability of porous media, p i is the pressure field of phase i, and g is the gravitational acceleration vector.
The pressure difference across the interface between the non-wetting phase (e.g., air) and the wetting phase (e.g., water) is called capillary pressure (p c ), and the pressure of the non-wetting phase is greater than the pressure of the wetting phase because of the curvature and interfacial tension of the interface (Brooks and Corey, 1964;[43, 44, 46]: where p s n is the pressure in the non-wetting phase, and p s w is the pressure in the wetting phase.
In the Brooks and Corey model, p c depends on phase saturation as shown in Eq. 4 Corey, 1964, 1966;[46]: where s w is the effective volume fraction of the wetting phase, which is expressed as: where s w is the volume fraction of the wetting phase, s rw is the residual volume fraction of the wetting phase, and s rn is the residual volume fraction of the non-wetting phase.
In the Brooks and Corey model, the relative permeability of the wetting phase ( k rs w ) and the non-wetting phase ( k rs n ) is expressed as Eqs. 6 and 7 Corey, 1964, 1966), respectively: where s n is the effective volume fraction of the nonwetting phase, which is expressed as: where s n is the volume fraction of the non-wetting phase.
In the "phtr" interface, there are two phases in the model: the water phase and the air phase. The "phtr" interface solves for the averaged volume fraction of water and air. The dependent variables are the volume fraction of water (s 1 ) and the volume fraction of air (s 2 ).
The sum of the two volume fractions follows the volume constraint equation: Equations 1 through 9 were implemented into the "phtr" interface of the COMSOL Multiphysics to simulate the water-air phase transport in overlay when exposed to weather variations.

Chloride transport
The transport of the chloride ions from high-concentration areas to low-concentration ones in the overlay obeys Fick's second law, which is expressed as Eq. 10 [9,19,23]: where C is chloride concentration, t is time, D is chloride diffusion coefficient, and x is a position variable.
In the "tds" interface of the COMSOL Multiphysics, Eq. 10 is written as Eqs. 11 and 12: where C i is the concentration of species i, t is time, J i is diffusive flux vector, R i is a reaction rate expression for species i (R i = 0 in this study), and D i is the diffusion coefficient of species i.

Model assumptions
To simplify the modeling process, the following assumptions were proposed for the two-dimensional axisymmetric model in this study: 1. Asphalt mastic, and the additives are homogeneous, isotropic, and incompressible material. 2. The overlay simulated in the model has been constructed and started to be used on May 1, 2012. 3. The subbase layer under the overlay is an impermeable asphalt layer. Thus, no water or chloride ions flows into the subbase layer. 4. The effect of gravity is neglected for water transport. 5. CaCl 2 is not deposited on the surface and pores of the overlay.

Model settings in COMSOL multiphysics
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% CaCl 2 -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 CaCl 2 -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 (P ec ) can be estimated using Eq. 13 [16,51]: where σ is interfacial tension of a fluid, θ is contact angle between the fluid and the capillary tube, and r c 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, r c is the average value of the pore's radii measured from scanning electron microscope (SEM) images. For the pores in asphalt mastic, the average r c was calculated with Eq. 14 [41]: where AV% is air void percentage or porosity in asphalt mastic, and D 75 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 P ec 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: where W i is the weight fraction of pores and W i = N i *A i /∑(N i *A i ), N i is the number of pores, and A i is the surface area of pores. N i and A i 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 (λ pmastic ) was calculated using Eq. 16 [33]: where k mastic 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/m 3 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 (s r ) 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 [m 2 ]' 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 (s 0,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:// weath er. 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]: where C(x,t) is the chloride concentration at a location x at a time t, C s is the chloride concentration at surface, C 0 is the initial chloride concentration at location x, erf(z) represents the error function: Thus, a steady-state chloride diffusion coefficient (D ss ) in asphalt mastic with a certain C s 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/m 3 ) of cement concrete with a watercement 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/m 3 , was introduced to represent chlorides potentially sourced from aggregates and binder, aimed to avoid the non-converged calculation Water permeability of epoxy layer on additives eps m 2 Fig. 4 a Relationship between D and C of asphalt pavement with CaCl 2 -zeolite/p-epoxy #16 and b goodness of fit between D from reference and calculated D with equation 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% CaCl 2 -zeolite/p-epoxy #16 [6,29,40,52]. C-related coefficient to D (A c ) is expressed as: where C s 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]: where D 0 is maximal diffusion coefficient, E a 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 (A T ) is expressed as: where T t 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 (19) adapted from [15], and it was implemented using the "Interpolation" function of COMSOL Multiphysics.
where T a is air temperature, R s 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 (22)  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 (C ave ) at the bottom of the simulation geometries over time was first determined with Eq. 23: where N n is the number of the non-cracking 150-mm section, C n is the chloride concentration at the bottom of the geometry of a non-cracking 150-mm section, N c is the number of the 1-crack 150-mm section, and C c is the chloride concentration at the bottom of the geometry of a 1-crack 150-mm section. C n and C c can be obtained from the simulation described below. Figure 6 shows the geometries used to determine C n and C c 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. the max crack depth of the overlay was 30 mm at 59%; 2. the max crack width was half of the max crack depth; 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 D ss 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 C s on the top and crack surfaces of the geometry was 560 mol/m 3 . The output was the relationship between the chloride concentration (e.g., C n and C c ) at the bottom of the geometry and time. Eleven points were selected to calculate C n and C c over time. The simulated C ave 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 C ave -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 (A tc ) is expressed as: The fatigue cracking-related coefficient to D (A fc ) was determined with the same method of A tc . 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 (24) A tc = e (2E−10)×t  Table S4 were obtained from the longterm pavement performance (LTPP) program of Federal Highway Administration, United States Department of Transportation.
It was assumed that the lane width of the 150 mmsection 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 C n and C c 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.
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, A fc is expressed as:  Combining D ss , Eq. 19, Eq. 21, Eq. 24, and Eq. 26, the final equation for the D of asphalt mastic (D mastic ) in the tds interface was expressed as: In addition, the threshold of moisture volume fraction (s 1 ) in concrete for chloride diffusion was 20% [35]. As such, D mastic was expressed as Eq. 27 when s 1 was equal or greater than 0.2, while D mastic was 0 when s 1 was less than 0.2, in which s 1 was the result of the "phtr" interface.
The parameters related to D mastic mentioned above are provided in Table 2.

Transport properties of additives
The D in the additives in the presence of water (D additive ) was estimated with the cooperation of the curve of (26) 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. D additive 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, D additive 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 m 2 /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 (C initial ) 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 D additive was 2E-11 m 2 /s and C initial was 1500 mol/m 3 . 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 D mastic , D additive was affected by s 1 from the result of the "phtr" interface. D additive was equal to 2E-11 m 2 /s when s 1 was equal or greater than 0.2, while D additive was 0 when s 1 was less than 0.2.

Initial values and boundary conditions
In reality, CaCl 2 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 s 1 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/m 3 ), 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/m 3 , while the initial value of the C in the asphalt mastic was assumed to be 0 mol/m 3 . 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

Materials
In this study, the same materials including asphalt binder and aggregates less than 4.75 mm, the same mixing design, and the same mixing and compaction procedure as in our previous study [53] were used to prepare asphalt mastic. The aggregates greater than 4.75 mm were not included in asphalt mastic. The air void of the asphalt mastic samples was 7% ± 0.5%. The same fabrication procedure and materials including CaCl 2 , zeolite, epoxy, halloysite nanoclay, and acetone as in our previous study [53] were used to prepare CaCl 2 -zeolite/p-epoxy #16. The calculated density of CaCl 2 -zeolite/p-epoxy #16 is 1.98 g/cm 3 .

Testing methods
To evaluate the diffusion behavior of chloride ions in asphalt mastic, a mastic diffusion test was adopted with a custom-designed apparatus as shown in Fig. 12. A mastic disk was sandwiched between two plastic cells, which were used to store the electrolyte solutions of CaCl 2 and Ca(NO 3 ) 2 , without any leakage. Mastic disk samples with a thickness of 12 mm ± 0.5 mm were cut from the prepared asphalt mastic samples and polished with sandpapers. The two cells were filled with DI water to saturate the mastic disk samples for 24 h before the mastic diffusion test, to minimize the effect of wick action (water transport) during ionic penetration [22]. Then, DI water was replaced by 0.28 mol/L (3 wt%) CaCl 2 electrolyte solution in the left cell and 0.28 mol/L Ca(NO 3 ) 2 electrolyte solution in the right cell. To avoid the evaporation of the water, the cells were covered by rubber caps. The variation of the chloride concentration over time in the right cell was periodically measured using a laboratory-prepared Ag/AgCl chloride sensor, which has been calibrated using standard NaCl solutions before each of measurements. The testing temperatures were 15 °C. The results of the test are shown in Figure S20. It is assumed that the D has been stable at the last point of the curve in Figure S20b, therefore D ss was equal to 1.14E-12 m 2 /s.
The inside-to-outside diffusion of chloride ions in CaCl 2 -zeolite/p-epoxy #16 was investigated by an additive immersion test, in which the manually compacted mixture of 133.6 g CaCl 2 -zeolite/p-epoxy #16 and 7.2 g asphalt (5.1 wt% binder content) was soaked in 1.5 L DI water at 15 °C as shown in Fig. 13. The mixing temperature was 145 °C, and the porosity of the mixture was 40%. The chloride concentration was monitored periodically by the same method used in the mastic diffusion test.
A falling head test [1], which is suitable for dense asphalt mixture samples with a low permeability, was used to evaluate k mastic . Figure 14 shows the schematic diagram of the setup of the falling head test. The diameter and thickness of mastic disk samples for the falling    where K mastic is the hydraulic conductivity of asphalt mastic, H 1 is the initial height of water in standpipe, H 2 is the final height of water in standpipe, A 1 is the cross-sectional area of sample, A 2 is the cross-sectional area of standpipe, L s is the length of sample, and t is the elapsed time of test. (29) where μ w is the dynamic viscosity of water at 25 °C, ρ w is the density of water, and g is the gravitational acceleration.
Water contact angle test was used to measure the values of θ of water to the epoxy layer and the asphalt mastic. An example of water contact angle tests on asphalt mastic is shown in Fig. 15. The results of the water contact angle test are shown in Table S8.
SEM (JEOL JXA8500F Field Emission Electron Microprobe) was adopted to measure the radii of the micropores in the epoxy coating layer of CaCl 2 -zeolite/p-epoxy. More details can be found in our previous study [53].

Results and discussion
Water transport Figures 16,17, and 18 depict the simulation results of the water transport model. Figure 16a-16f display the representative images of the evolution of the s 1 in the overlay when exposed to varying precipitation and temperature. At the very beginning (Fig. 16a, t = 0 h), the s 1 regardless of the location was in the model was equal to 0 because there was no rainfall or snowfall on the overlay. At 72 h (Fig. 16b), the s 1 increased to different levels, depending on the location in the model. Two rainfalls between 0 and 72 h accounted for the increase in the s 1 , which could be observed in Fig. 17. Figure 17a illustrates the temporal evolution of the s 1

Fig. 15
Representative water contact angle test on asphalt mastic  Figure 17b reveals that the first rainfall appeared between 56 and 60 h and the second one rained between 63 to 65 h. In contrast, the drying process between 60 and 63 h and after 65 h was responsible for the different values of the s 1 at the different locations. Water could only evaporate through the top surface, which explains why the s 1 near the top surface had lower values than that around the middle and the bottom. At 96 h (Fig. 16c), the s 1 changed slightly due to the corporation of the third 1-hour rainfall (as shown in Fig. 17b) and the drying process after 90 h. Between 96 and 480 h, there was no rainfall. Therefore, the s 1 gradually decreased to low values over time, as shown in Fig. 16d. Figure 16e reveals that all the s 1 in the overlay reached to 1 at 840 h, resulting from sufficient amount of rains falling on the surface before 840 h. The last status of the overlay about s 1 for this study was shown in Fig. 16f (44,280 h, about 5.05 years). Figure 18 illustrates the temporal evolution of the s 1 at the center of the geometry (in asphalt mastic), which was consistent with the results shown in Figs. 16 and 17.
It should be noted that pavement temperature is a variable in several equations of the model's inputs, including the entry capillary pressure of water in pores of asphalt mastic and functional additive, the dynamic viscosity of water and air, and the density of air. Pavement temperature affects the water fraction through influencing these parameters. Chloride Transport Figure 19 displays the representative images that illustrate the temporal evolution of the C in the overlay when exposed to wet-dry cycles, varying pavement temperature, pavement thermal cracking, and pavement fatigue cracking. The wet-dry cycles in the model were expressed by s 1 (> = 0.2 for a wet condition and < 0.2 for a dry condition), which resulted from the prediction of the water transport model. The C in asphalt mastic was 0 mol/m 3 from 0 to 48 h as shown in Fig. 19a-c, because the overlay was in a dry condition (s 1 < 0.2) as shown in Fig. 17b. After the first rainfall, the s 1 in asphalt mastic was greater than 0.2, thus the C gradually increased over time (as shown in Fig. 19d-f ) until the s 1 decreased to be lower than 0.2 due to water evaporation (when there was no rainfall on the pavement). Figure 19f-h indicate that the C in asphalt mastic remained unchanged from 144 to 240 h. The reason was that the D mastic was equal to 0 when the s 1 was lower than 0.2 during that period as shown in Fig. 18b. At 4272 h (Fig. 19i), the C in asphalt mastic greatly increased while the C in functional additives significantly decreased, relative to those before 240 h. Chloride ions in the overlay moved gradually from the locations with high concentrations to the locations with low concentrations. The last status of the C in asphalt mastic and functional additive for this study was shown in Fig. 19j (44,280 h, about 5.05 years). In this figure, the C in asphalt mastic near the top surface (the blue area) was visibly lower than that at the middle and bottom, which is due to dilution and removal of chloride ions from the top surface by precipitation of rain or snow. The C at deeper areas has higher values, mainly because the chloride ions in the additives at deeper areas need to travel longer distances to reach the top surface and their transport takes more time than the time needed by the dilution.
In addition to moisture availability (due to precipitations), pavement temperature, pavement thermal cracking, and pavement fatigue cracking were three coefficients in the equation of the diffusion coefficient of asphalt mastic, all of which influenced the change of the C in asphalt mastic and functional additive accordingly.
Although Fig. 19 qualitatively illustrates the overall variation of the C in the overlay, it is desirable to quantitatively monitor the continuous evolution of the C in functional additive and asphalt mastic and on the top surface, respectively. Therefore, Figs. 20, 21, and 22 were plotted accordingly. Figure 20 contrasts the temporal evolution of the C in three functional additives at three different locations: top, middle, and bottom. They all shared a same trend and gradually decreased over time. In some periods the curve of the C exhibited a stable plateau phase due to that the Ds of asphalt mastic and functional additive were 0 m 2 /s because the s 1 was equal to 0. The slight differences among the three Cs resulted from the surrounding environment of the three additives, such as the quantity of functional additives and aggregates nearby, the shape of Fig. 19 Representative images of evolution of C in APSSA when exposed to wet-dry cycles and variations of temperature, pavement thermal cracking, and pavement fatigue cracking. From a to j, t = 0, 24,72,96,144,192,240,4272 adjacent aggregates, and the cumulative time when s 1 was equal or greater than 0.2. Figure 21 illustrates the temporal evolution of the C on the top surface of the overlay. Five points on the top surface were selected, and the averaged C was plotted. It can be seen that the C on the top surface was reduced following the pattern of precipitations shown in Fig. 17. The amount of the reduction in C depended on the type of precipitation; the reduction in C was not exact 95% and 50% for rainfalls and snowfalls respectively. This is because that the setting in COMSOL Multiphysics for the chloride removal at top surface was the boundary condition 'Flux' , instead of directly removing the available C. Arguably, the plot shown in Fig. 21 better represents the reality.
The effective anti-icing concentration of Cl − in water for different temperatures is provided in Table 3, which was converted based on the phase diagram for the water-CaCl 2 system as shown in Figure S21. One can assume that there is a linear relationship between temperature and the CaCl 2 concentration (in wt%) when the temperature is higher than -10 °C. The equation is expressed as: where C CaCl2 is the effective mass percentage of CaCl 2 for anti-icing at T and T is temperature. The temperature of asphalt pavement near surface is usually 1 to 2 °C higher than air temperature during winter [34], and snow events generally occur when air temperature is lower than 0 °C in winter. As such, this study assumes that snow events occurred when the temperature of the overlay dropped below 2 °C, and the historical hourly pavement temperature data (near Pullman, WA) below 2 °C were thus used to calculate the statistical parameters (e.g., percentiles mentioned later).
Considering the compressed or vibrated pavement surface due to the effects of loading and friction from moving vehicles [54] and the average size of functional additive being 2 mm, this study assumes that the C in the overlay at a depth of 4 mm approximately determines whether the anti-icing function is effective or not. The averaged C of five points in asphalt mastic at the depth of 4 mm was plotted as shown in Fig. 22. At the beginning, the C increased rapidly and reached to 1,118 mol/m 3 within two months. The fast diffusion resulted from the big difference of the C between asphalt mastic and functional additives. Subsequently, the C gradually decreased over time, due to chloride removal on the top surface. According to the plot in Fig. 22 and the data in Table 3, the anti-icing function of the overlay was fully effective in 17,520 h (2 years) and 43,800 h (5 years) for the minimum pavement temperature above -3.4 °C and -2.4 °C, corresponding to its lowest C of 808 mol/m 3 and 569 mol/m 3 , respectively. These two pavement temperatures translate to 97.4-percentile and 96.3-percentile of historical hourly pavement temperature near Pullman, Washington, respectively. In other words, this specific 16-mm thick functional overlay could prevent the surface from icing, for 97.4% of the snow events within 2 years and for 96.3% of the snow events within 5 years. The simulation results agree well with the estimated anti-icing life using Eq. 3.2 in our previous study [53]. It should be noted that these effective anti-icing service lives are conservatively calculated, because beyond them there are still residual friction benefits of this functional overlay, because part of the overlay surface would still be able to keep the moisture in a liquid instead of frozen state.

Sensitivity analysis
Sensitivity analysis was conducted to investigate how the C of the overlay was affected by variations in the overlay thickness (h) and D mastic , because the simulation  Figure 23a depicts the variation of the C over time at 4 mm below the top surface as a function of the h and D mastic settings, over the simulated 43,800 h (5 years). In the first 1,500 h, as shown in Fig. 23b, the C decreased with reducing h or D mastic , and a combined action of reducing both h and D mastic further decreased the C. In contrast, the C increased with increasing h or D mastic , and the C further increased when h and D mastic increased simultaneously. For a long-term evaluation (as shown in Fig. 23a), interestingly, the reduction of D mastic or both D mastic and h resulted in slight increase in the C on the top surface of the overlay, likely due to a smaller amount of chloride ions being removed by precipitations. The other six changes of the parameters slightly decreased the C, indicative of multiple mechanisms at work. Table 4 shows the predicted anti-icing life of the functional overlay for the different settings, assuming the goal is to prevent the surface from icing for 96.3% of the snow events. Under the original setting (h = 16 mm, D mastic = Eq. 27), the anti-icing function of the overlay is predicted to be effective for the overlay temperature above -2.4 °C, and the corresponding Cl − concentration was 557 mol/m 3 , as shown in Table 3. Compared with the original setting, the anti-icing life of the overlay increased by 2% when D mastic decreased by 20% and increased by 1% when both D mastic and h decreased by 20%. On the other hand, it decreased by 6%, 9%, 9%, 10%, 1%, and 11% when D mastic increased by 20%, h decreased by 20%, h increased

Indirect validation of the predictive model
This study is limited to exploratory investigation. The finite element method-based model was built by following first principles of water and chloride transport, incorporated experimental results from laboratory, and relying on case study data and assumptions where necessary. For the prediction of anti-icing life of the overlay with the functional additive, direct validation of this model remains challenging due the lack of directly usable data in the published domain and the lack of resources to conduct long-term laboratory or field investigation for longevity validation. Nonetheless, this section describes encouraging results that indirectly validate the predictive model established thus far.   Figure 24 displays the temporal evolution of water fraction on the top surface vs. at center point of the overlay, and temporal evolution of chloride concentration at center point of the overlay vs. in an additive near center point, up to 1500 h. The variations of water fraction on the top surface of the overlay closely matches the pattern of how precipitation occurred on an hourly basis. It can be clearly seen that the water fraction at the center point rapidly increases during precipitations, whereas it gradually decreases due to natural water evaporations when there is no precipitation. Moreover, the chloride concentration at the center point increases over time when the water fraction is equal to or greater than 0.2, whereas the chloride concentration remains unchanged when the water fraction is less than 0.2. In contrast, the chloride concentration in the additive decreases over time when the water fraction is equal to or greater than 0.2, whereas it remains unchanged when the water fraction is less than 0.2. These predicted trends for both the water fraction and the chloride concentration at different points in the overlay are reasonable and consistent with the reality.
In our previous study [53], the results of the fog-freezing tests indicated the friction coefficient of the asphalt mixture with 5.1 wt% additive #16 decreased from 0.5 at -3.9 °C to 0.4 at -9.4 °C. This can be explained by the temporal evolution of the chloride concentration at a depth of 4 mm below the surface of the overlay. Figure 22 illustrates the predicted profile of the chloride concentration in the functional overlay at 1248 h. The model results indicate that the lowest temperature for an effective anti-icing of this overlay was about -4.6 °C (corresponding to its highest chloride concentration of 1,118 mol/m 3 on the plot, see Fig. 23a). At this temperature or above, all the water remained a liquid state and would not turn into ice. Thus, the anti-icing function was fully effective at -3.9 °C, which explains the high friction coefficients from the fog-freezing experiment at that temperature. It is apparent that 1,118 mol/m 3 would not be sufficient to keep all the water (at the overlay surface) at a liquid state at -9.4 °C; instead, a minimum C of 2,400 mol/m 3 is needed, by the theoretical evaluation based on the eutectic curve (see Figure  S21). The simulation result at 1248 h revealed the chloride concentration at some locations at the 4-mm depth was higher than 2,400 mol/m 3 , as shown in the areas marked with purple ellipses in Fig. 25. These localized high chloride-concentration locations help explain the experimental result at -9.4 °C. Because not all the water on the sample surface was prevented from icing, the friction coefficient reduced to 0.4 at -9.4 °C from 0.5 at -3.9 °C.

Conclusions
This study employed a Finite Element Method based software to predict the anti-icing life of a thin overlay of asphalt pavement with salt-storage additive, by simulating chloride transport in the overlay when exposed to the changes in weather conditions in terms of precipitation and temperature and pavement conditions in terms of thermal cracking and fatigue cracking. The FEM based simulation included water transport modeled with the "phtr" interface and chloride transport modeled with the "tds" interface of the software. The output of the water transport model was water fraction, which was used as one input parameter of the chloride transport model. The water fraction plays a crucial role in controlling the chloride diffusion coefficients of asphalt mastic and functional additive.
The results of the water transport model show the temporal and spatial evolutions of the water fraction in the overlay as a function of (historical) precipitation and pavement temperature. The water fraction increased rapidly during precipitations, whereas it gradually decreased due to water evaporation through the surface of the overlay when there was no precipitation. It was easier to lose water at the locations near the surface of the overlay than that around the middle and bottom. Pavement temperature was a variable affecting many parameters of the model, thus affecting the water fraction. These parameters included entry capillary pressure of water in pores of asphalt mastic and functional additive, dynamic viscosity of water and air, and density of air.
For the chloride transport model, the chloride concentration in asphalt mastic highly depended on the value of the water fraction, i.e., moisture availability limited the Cl − diffusion coefficients of asphalt mastic and functional additive. The two Cl − diffusion coefficients were equal to 0 when the water fraction was lower than 0.2. As the water fraction was equal to or higher than 0.2, the Cl − diffusion coefficient of asphalt mastic was defined as a function of the surrounding chloride concentration, pavement temperature, pavement thermal cracking, and pavement fatigue cracking, and the Cl − diffusion coefficient of functional additive remained constant. Chloride ions migrated from high-concentration locations to lowconcentration ones when the water fraction was equal to or higher than 0.2, until a balanced concentration profile was obtained in the system. The chloride concentration in asphalt mastic increased over time and the chloride concentration in functional additive decreased over time when the water fraction was equal to or higher than 0.2, otherwise it remained unchanged. Due to chloride removal on the top surface, the chloride concentration near the top surface was lower than that around the middle and bottom of the overlay. Although the model considered the effect of pavement thermal cracking and fatigue cracking on increasing the Cl− diffusion coefficient of asphalt mastic, other pavement distresses (moisture damage, rutting, longitudinal cracking, UV aging, etc.) and the interactions among these distresses were not included. Additional research is
Additional file 1: Figure S1. Relationship between σ of water and pavement temperatures. Figure S2. Particle-size distribution curve of aggregates used for asphalt mixture. Figure S3. Relationship between μ of water and pavement temperatures. Figure S4. Relationship between μ of air and pavement temperatures. Figure S5. Relationship between ρ of air and pavement temperatures. Figure S6. a) Simulated and calculated C-time curves and b) goodness of fit at 0% cracking level for thermal cracking. Figure S7. a) Simulated and calculated C-time curves and b) goodness of fit at 27% cracking  Figure S8. a) Simulated and calculated C-time curves and b) goodness of fit at 38% cracking level for thermal cracking. Figure S9. a) Simulated and calculated C-time curves and b) goodness of fit at 59% cracking level for thermal cracking. Figure S10. ESALs-fatigue crack rate relationship of asphalt pavement. Figure S11. EASLs-crack width relationship of asphalt pavement. Figure S12. Loading number-crack width relationship of asphalt pavement. Figure S13. Load-single LEF relationship of asphalt pavement. Figure S14. Crack width-crack depth relationship of asphalt pavement. Figure S15. a) Simulated and calculated C-time curves and b) goodness of fit at 0 year for fatigue cracking. Figure S16. a) Simulated and calculated C-time curves and b) goodness of fit at 8 years for fatigue cracking. Figure S17. a) Simulated and calculated C-time curves and b) goodness of fit at 9 years for fatigue cracking. Figure S18. a) Simulated and calculated C-time curves and b) goodness of fit at 10 years for fatigue cracking. Figure S19. a) C-time curves from simulation and experiment and b) goodness of fit to determine D of additives (datapoints before 71 h were removed in Figure S19b due to unstable diffusion). Figure S20. a) Cl concentration-time relationship and b) Cl diffusion coefficient-time relationship from mastic diffusion test. Figure S21. Phase diagram for water-CaCl2 system (adopted from Phase diagrams for binary salt solutions | Phasediagram). Table S1. a) Area fractions of asphalt mastic, aggregates, and additives in each of samples. Table S2. An example of calculating λ p-epoxy for Figure 2. Table S3. Estimated λ p-epoxy based on three SEM images of the surface of CaCl 2 -zeolite/p-epoxy #16.