Finite Element Method-Based Skid Resistance Simulation Using In-Situ 3D Pavement Surface Texture and Friction Data

Skid resistance is an important surface characteristic that influences roadway safety. Various studies have been performed to understand the interaction between pavement and tires through numerical simulation for skid resistance prediction. However, the friction parameters required for simulation inputs are generally determined by objective assumptions. This paper develops a finite element method (FEM)-based skid resistance simulation framework using in-situ 3D pavement surface texture and skid resistance data. A 3D areal pavement model is reconstructed from high resolution asphalt pavement surface texture data. The exponential decay friction model is implemented in the simulation and the interface friction parameters required for the simulation are determined using the binary search back-calculation approach based on a trial process with the desired level of differences between simulated and observed skid numbers. To understand the influence of texture characteristics on interface friction parameters, the high-resolution 3D texture data is separated into macro- and micro-scales through Butterworth filtering and various areal texture indicators are calculated at both levels. Principal component analysis (PCA) regression analysis is conducted to quantify the relationship between various texture characteristics and the interface friction parameters. The results from this study can be used to better prepare the inputs of friction parameters for FEM simulation.


Introduction
The risk of traffic accidents may rise significantly when the skid resistance is lower than a certain threshold [1]. The skid resistance is measured as a resistive drag force, generally using the locked-wheel, dynamic friction tester or grip tester with standard testing tires or rubber sliders [2]. The pavement interface friction is affected by many factors, such as vehicle factors (load, speed, slip ratio), rubber properties, asphalt pavement factors (aggregate shape, roughness, micro-and macro-texture), and weather conditions (temperature and contamination) [3][4][5].
Many researchers have contributed to monitoring and predicting pavement friction in the past decades [6]. Macro-and micro-pavement textures have been found to contribute significantly to surface friction and various relationships have been developed [3,[7][8][9]. With advances in noncontact three-dimensional (3D) measurement technologies and developments in high performance computers, wavelet analysis, the Hilbert-Huang transform, fractal analysis, power spectra density, and Persson's model have been used to characterize pavement macrotexture attributes and correlate them with

Objective
The objective of this paper is to develop a rubber pavement interaction simulation framework to determine the interface friction of pavement surfaces using in-situ 3D pavement surface texture and skid resistance value. As illustrated in Figure 1, the framework is composed of the following components: • Collecting in-situ high-resolution 3D pavement surface texture and skid resistance measured by the dynamic friction tester (DFT); • Characterizing pavement texture attributes at both macro-and micro-scale with 3D areal texture parameters; the Butterworth filter is applied to separate the texture data into two scales; • Proposing a re-construction method to establish 3D areal pavement surfaces as inputs for the FEM numerical simulation below; • Implementing the FEM-based rubber pavement interaction model to simulate the DFT measurements; • Computing the rubber pavement interface friction using the DFT data sets according to the binary search back-calculation method; • Determining the relationships between rubber pavement interface friction and the 3D areal pavement texture parameters using principal component analysis PCA regression.

Field Data Collection
The field-testing bed of this study is the long-term pavement performance (LTPP) Specific Pavement Study 10 (SPS-10) sites, which were constructed by the Oklahoma Department of Transportation (DOT) on Highway 66 in Yukon in November 2015. The experimental matrix, to evaluate the short-and long-term performance of the asphalt mixtures, includes one hot mix asphalt (HMA) and two warm-mix asphalt (WMA) experimental treatment sections, where a foaming process and a chemical additive with 10% to 25% reclaimed asphalt pavement (RAP) and reclaimed asphalt shingle (RAS) content are required per LTPP [39]. Under the SPS-10 experiment treatment, six LTPP SPS-10 experimental treatment sections and six control sections were constructed in Oklahoma. The site location and the corresponding length for each section are shown in Figure 2. The detailed mixture design for the sites is presented by Table 1. High resolution pavement 3D surface texture and DFT skid resistance data sets were collected in the field on the same day in January 2017. The average ambient temperature during the testing is 16 °C.

Field Data Collection
The field-testing bed of this study is the long-term pavement performance (LTPP) Specific Pavement Study 10 (SPS-10) sites, which were constructed by the Oklahoma Department of Transportation (DOT) on Highway 66 in Yukon in November 2015. The experimental matrix, to evaluate the short-and long-term performance of the asphalt mixtures, includes one hot mix asphalt (HMA) and two warm-mix asphalt (WMA) experimental treatment sections, where a foaming process and a chemical additive with 10% to 25% reclaimed asphalt pavement (RAP) and reclaimed asphalt shingle (RAS) content are required per LTPP [39]. Under the SPS-10 experiment treatment, six LTPP SPS-10 experimental treatment sections and six control sections were constructed in Oklahoma. The site location and the corresponding length for each section are shown in Figure 2. The detailed mixture design for the sites is presented by Table 1. High resolution pavement 3D surface texture and DFT skid resistance data sets were collected in the field on the same day in January 2017. The average ambient temperature during the testing is 16 • C.

Field Data Collection
The field-testing bed of this study is the long-term pavement performance (LTPP) Specific Pavement Study 10 (SPS-10) sites, which were constructed by the Oklahoma Department of Transportation (DOT) on Highway 66 in Yukon in November 2015. The experimental matrix, to evaluate the short-and long-term performance of the asphalt mixtures, includes one hot mix asphalt (HMA) and two warm-mix asphalt (WMA) experimental treatment sections, where a foaming process and a chemical additive with 10% to 25% reclaimed asphalt pavement (RAP) and reclaimed asphalt shingle (RAS) content are required per LTPP [39]. Under the SPS-10 experiment treatment, six LTPP SPS-10 experimental treatment sections and six control sections were constructed in Oklahoma. The site location and the corresponding length for each section are shown in Figure 2. The detailed mixture design for the sites is presented by Table 1. High resolution pavement 3D surface texture and DFT skid resistance data sets were collected in the field on the same day in January 2017. The average ambient temperature during the testing is 16 °C.
Note: E1-E6 the experimental treatment sections; C1-C6 the control sections.    In particular, two instruments are used in the field for texture and skid resistance data collection. 3D pavement surface texture data was obtained utilizing the LS-40 porTable 3D surface analyzer, as shown in Figure 3. This analyzer scans an area of 101.6 mm by 114.3 mm with the height resolution (z) of 0.01 mm and lateral resolution (x,y) of 0.05 mm [40]. The high-resolution 3D texture data acquired from LS-40 includes both macro-and micro-level texture information of the scanned surfaces. Each scan has 2048 by 2448 cloud points. In particular, two instruments are used in the field for texture and skid resistance data collection. 3D pavement surface texture data was obtained utilizing the LS-40 porTable 3D surface analyzer, as shown in Figure 3. This analyzer scans an area of 101.6 mm by 114.3 mm with the height resolution (z) of 0.01 mm and lateral resolution (x,y) of 0.05 mm [40]. The high-resolution 3D texture data acquired from LS-40 includes both macro-and micro-level texture information of the scanned surfaces. Each scan has 2,048 by 2,448 cloud points.  Skid resistance data was collected using DFT, which measures pavement surface frictional properties at various speeds [41]. DFT is widely used for skid resistance because it is repeatable and reproducible with controlled operating procedures or ambient factors [42]. It consists of a horizontal spinning disk mounted with three spring-loaded rubber sliders, as shown in Figure 4. Water spray in front of the rubber sliders to form 1 mm water film thickness and a constant vertical load is applied on the slider, while the disk spins on the test surface. The torque signal is monitored constantly while the velocity of the spinning disk decreases due to the friction between the rubber sliders and the test surface, thus the pavement surface skid resistance data is derived.
Within each LTPP SPS-10 section, three pairs of LS-40 3D surface texture data and DFT skid resistance data were obtained at 100-ft intervals, starting at the beginning of each section. On the mainline control section after each experimental treatment section, an additional three pairs of pavement texture and skid resistance measurements were conducted at 300-ft intervals. Therefore, 36 pairs of 3D pavement texture and skid resistance data measurements were obtained for each data collection.

Interface Friction Model
Interface friction occurs at the contact surface between the pavement and the tire, resulting from adhesion and hysteresis. The adhesion is related to interface shear strength while the hysteresis is the result of damping losses and energy dissipation of the rubber arising from the pavement surface asperities [32].  Skid resistance data was collected using DFT, which measures pavement surface frictional properties at various speeds [41]. DFT is widely used for skid resistance because it is repeatable and reproducible with controlled operating procedures or ambient factors [42]. It consists of a horizontal spinning disk mounted with three spring-loaded rubber sliders, as shown in Figure 4. Water spray in front of the rubber sliders to form 1 mm water film thickness and a constant vertical load is applied on the slider, while the disk spins on the test surface. The torque signal is monitored constantly while the velocity of the spinning disk decreases due to the friction between the rubber sliders and the test surface, thus the pavement surface skid resistance data is derived.
Within each LTPP SPS-10 section, three pairs of LS-40 3D surface texture data and DFT skid resistance data were obtained at 100-ft intervals, starting at the beginning of each section. On the mainline control section after each experimental treatment section, an additional three pairs of pavement texture and skid resistance measurements were conducted at 300-ft intervals. Therefore, 36 pairs of 3D pavement texture and skid resistance data measurements were obtained for each data collection. In particular, two instruments are used in the field for texture and skid resistance data collection. 3D pavement surface texture data was obtained utilizing the LS-40 porTable 3D surface analyzer, as shown in Figure 3. This analyzer scans an area of 101.6 mm by 114.3 mm with the height resolution (z) of 0.01 mm and lateral resolution (x,y) of 0.05 mm [40]. The high-resolution 3D texture data acquired from LS-40 includes both macro-and micro-level texture information of the scanned surfaces. Each scan has 2,048 by 2,448 cloud points.  Skid resistance data was collected using DFT, which measures pavement surface frictional properties at various speeds [41]. DFT is widely used for skid resistance because it is repeatable and reproducible with controlled operating procedures or ambient factors [42]. It consists of a horizontal spinning disk mounted with three spring-loaded rubber sliders, as shown in Figure 4. Water spray in front of the rubber sliders to form 1 mm water film thickness and a constant vertical load is applied on the slider, while the disk spins on the test surface. The torque signal is monitored constantly while the velocity of the spinning disk decreases due to the friction between the rubber sliders and the test surface, thus the pavement surface skid resistance data is derived.
Within each LTPP SPS-10 section, three pairs of LS-40 3D surface texture data and DFT skid resistance data were obtained at 100-ft intervals, starting at the beginning of each section. On the mainline control section after each experimental treatment section, an additional three pairs of pavement texture and skid resistance measurements were conducted at 300-ft intervals. Therefore, 36 pairs of 3D pavement texture and skid resistance data measurements were obtained for each data collection.

Interface Friction Model
Interface friction occurs at the contact surface between the pavement and the tire, resulting from adhesion and hysteresis. The adhesion is related to interface shear strength while the hysteresis is the result of damping losses and energy dissipation of the rubber arising from the pavement surface asperities [32].

Interface Friction Model
Interface friction occurs at the contact surface between the pavement and the tire, resulting from adhesion and hysteresis. The adhesion is related to interface shear strength while the hysteresis is the result of damping losses and energy dissipation of the rubber arising from the pavement surface asperities [32].
Many studies have concluded that the skid resistance performance mainly depends on the rubber pavement interface friction [2,3], as illustrated in Figure 5. To describe the interface friction property, the exponential decay friction model proposed by Oden and Martins [34] is generally used. One example DFT testing is illustrated in Figure 6, which has demonstrated a similar trend as that of the exponential decay friction model. The exponential decay function is shown in Equation (1): where µ s is the static friction coefficient, µ k is the kinetic friction coefficient, α is a user-defined decay coefficient, and s is the sliding velocity. In this study, the rubber pavement interaction is lubricated with water films during the simulation [43]. The influence of water on the interface friction is reflected in the three friction input parameters in the exponential decay function: µ s the static friction coefficient, µ k the kinetic friction coefficient, and α the decay coefficient. Many studies have concluded that the skid resistance performance mainly depends on the rubber pavement interface friction [2,3], as illustrated in Figure 5. To describe the interface friction property, the exponential decay friction model proposed by Oden and Martins [34] is generally used. One example DFT testing is illustrated in Figure 6, which has demonstrated a similar trend as that of the exponential decay friction model. The exponential decay function is shown in Equation (1): where is the static friction coefficient, is the kinetic friction coefficient, is a user-defined decay coefficient, and is the sliding velocity. In this study, the rubber pavement interaction is lubricated with water films during the simulation [43]. The influence of water on the interface friction is reflected in the three friction input parameters in the exponential decay function: the static friction coefficient, the kinetic friction coefficient, and α the decay coefficient.

Reconstruction of 3D Pavement Texture Surface
To accurately illustrate the rubber pavement contact mechanism, many pavement surface models have been established based on sine patterns [44], X-ray tomography [36], simplified porous pavement surface [45], hemispheric roughness surface [46] and other forms [47][48][49][50], to reveal the transient dynamic performance of rubber when the tire traversed over a pavement segment [44].
In this paper, the areal pavement surface model is reconstructed based on the high-resolution 3D LS-40 texture data sets. Figure 7 shows the reconstruction procedure of the 3D texture surface for FEM simulation mesh using field data. Pavement surface texture images are obtained using LS-40 (HyMIT Measurement Instrument Technology, Austin, TX, USA) from the field and saved as 2048 × 2448 16-bit range data. Speckle noises can exist in the LS-40 pavement texture range data. Subsequently, noises are eliminated by applying the limiting filter algorithm [51] as following: for the cloud points in each pavement profile frame, the first quartile, median, third quartile are firstly Many studies have concluded that the skid resistance performance mainly depends on the rubber pavement interface friction [2,3], as illustrated in Figure 5. To describe the interface friction property, the exponential decay friction model proposed by Oden and Martins [34] is generally used. One example DFT testing is illustrated in Figure 6, which has demonstrated a similar trend as that of the exponential decay friction model. The exponential decay function is shown in Equation (1): where is the static friction coefficient, is the kinetic friction coefficient, is a user-defined decay coefficient, and is the sliding velocity. In this study, the rubber pavement interaction is lubricated with water films during the simulation [43]. The influence of water on the interface friction is reflected in the three friction input parameters in the exponential decay function: the static friction coefficient, the kinetic friction coefficient, and α the decay coefficient.

Reconstruction of 3D Pavement Texture Surface
To accurately illustrate the rubber pavement contact mechanism, many pavement surface models have been established based on sine patterns [44], X-ray tomography [36], simplified porous pavement surface [45], hemispheric roughness surface [46] and other forms [47][48][49][50], to reveal the transient dynamic performance of rubber when the tire traversed over a pavement segment [44].
In this paper, the areal pavement surface model is reconstructed based on the high-resolution 3D LS-40 texture data sets. Figure 7 shows the reconstruction procedure of the 3D texture surface for FEM simulation mesh using field data. Pavement surface texture images are obtained using LS-40 (HyMIT Measurement Instrument Technology, Austin, TX, USA) from the field and saved as 2048 × 2448 16-bit range data. Speckle noises can exist in the LS-40 pavement texture range data. Subsequently, noises are eliminated by applying the limiting filter algorithm [51] as following: for the cloud points in each pavement profile frame, the first quartile, median, third quartile are firstly calculated. If the values of the range data are greater than 1.5 times of the third quartile value, or smaller than 1.5 times of the first quartile, they are treated as noises, whose values are replaced with the third quartile or the first quartile.

Reconstruction of 3D Pavement Texture Surface
To accurately illustrate the rubber pavement contact mechanism, many pavement surface models have been established based on sine patterns [44], X-ray tomography [36], simplified porous pavement surface [45], hemispheric roughness surface [46] and other forms [47][48][49][50], to reveal the transient dynamic performance of rubber when the tire traversed over a pavement segment [44].
In this paper, the areal pavement surface model is reconstructed based on the high-resolution 3D LS-40 texture data sets. Figure 7 shows the reconstruction procedure of the 3D texture surface for FEM simulation mesh using field data. Pavement surface texture images are obtained using LS-40 (HyMIT Measurement Instrument Technology, Austin, TX, USA) from the field and saved as 2048 × 2448 16-bit range data. Speckle noises can exist in the LS-40 pavement texture range data. Subsequently, noises are eliminated by applying the limiting filter algorithm [51] as following: for the cloud points in each pavement profile frame, the first quartile, median, third quartile are firstly calculated. If the values of the range data are greater than 1.5 times of the third quartile value, or smaller than 1.5 times of the first quartile, they are treated as noises, whose values are replaced with the third quartile or the first quartile.
Five categories of 3D areal texture parameters are calculated at the macro-and micro-texture scales and for the raw images before the separation: height parameters, spatial parameters, hybrid parameters, volume parameters, and feature parameters [40]. All the texture parameters are processed via the MountainsMap ® software package (Version 7.3, Digital Surf, Besançon, Bourgogne-Franche-Comté, France). The relationships between the interface friction and 3D areal texture parameters for the macro-and micro-texture data sets are discussed later.

FEM based Skid Resistance Simulation
In this study it is assumed that the deformations of the pavement surface are negligible as compared to that of the rubber. The 3D areal pavement surface is re-constructed as a rigid body to increase computational efficiency. The size of the pavement model is 114.3 mm × 40 mm × 8 mm.
Since sharp angles can lead to computing error and convergence problems [54], the pavement element size is set as 0.5 mm. In addition, surface smoothing techniques are applied before meshing pavement surfaces. Rubber sliding on pavements is a transient dynamic behavior. During the DFT testing, the rubber block is pressing and sliding against the pavement surface. Large deformation could occur during the process, which could result in convergence problems in simulation. Therefore, the 3D linear eight-node brick element (C3D8R) [55] is selected to model the rubber slider of DFT. The size of the rubber block is 20 mm × 16 mm × 6 mm [41].
A rubber slider is a near incompressible and hyper-elastic material with viscoelasticity [56]. It is synthetic rubber as specified in the ASTM E501 specification [41,57]. The DFT rubber slider has the The influence of surface texture on friction occurs at both the macro-and micro-scale [14]. In this study, the acquired LS-40 texture data is separated into macro-and micro-texture data by using the Butterworth filter [53]. All the frequencies between 0.0008 and 0.08 cycles/m (wavelengths from 0.5 to 50 mm) are passed to isolate only the effect of macro-texture, while all the frequencies less than 0.08 cycles/m (wavelengths lower than 0.5 mm) are saved to represent the micro-texture information. Five categories of 3D areal texture parameters are calculated at the macro-and micro-texture scales and for the raw images before the separation: height parameters, spatial parameters, hybrid parameters, volume parameters, and feature parameters [40]. All the texture parameters are processed via the MountainsMap ® software package (Version 7.3, Digital Surf, Besançon, Bourgogne-Franche-Comté, France). The relationships between the interface friction and 3D areal texture parameters for the macroand micro-texture data sets are discussed later.

FEM Based Skid Resistance Simulation
In this study it is assumed that the deformations of the pavement surface are negligible as compared to that of the rubber. The 3D areal pavement surface is re-constructed as a rigid body to increase computational efficiency. The size of the pavement model is 114.3 mm × 40 mm × 8 mm. Since sharp angles can lead to computing error and convergence problems [54], the pavement element size is set as 0.5 mm. In addition, surface smoothing techniques are applied before meshing pavement surfaces. Rubber sliding on pavements is a transient dynamic behavior. During the DFT testing, the rubber block is pressing and sliding against the pavement surface. Large deformation could occur during the process, which could result in convergence problems in simulation. Therefore, the 3D linear eight-node brick element (C3D8R) [55] is selected to model the rubber slider of DFT. The size of the rubber block is 20 mm × 16 mm × 6 mm [41].
A rubber slider is a near incompressible and hyper-elastic material with viscoelasticity [56]. It is synthetic rubber as specified in the ASTM E501 specification [41,57]. The DFT rubber slider has the same compounding requirements as the standard tires used for common pavement skid-resistance testing devices, such as locked-wheel trailer and grip tester [41,[57][58][59]. However, rubber manufacturers usually do not publish material's mechanical property information. Many previous research studies have simulated the tire-pavement interaction by setting the rubber as a linear elastic material with a Poisson ratio around 0.5 [30,32,45,60]. It is also pointed out that the linear elastic material model still has the potential to demonstrate the rubber's mechanical behavior by [61]. Hence, in this research, the linear elastic material properties are adopted from the papers by Ong et al. and Zhang et al. [45,60], in which the tire rubber material property is as same as the DFT rubber slider's. The elastic modulus of the rubber slider is 100 MPa and the Poisson's ratio and density of the rubber are set as 0.45 and 1200 kg/m 3 , and rubber's viscoelastic parameter's is adopted from the papers by Zegard et al. [61], respectively. The applied load on the slider is 11.8 N, and the sliding speeds for simulation are 20, 40 and 60 km/h. Since all the data collection was completed with very similar temperature conditions, the properties of the rubber material are considered to be constant. The rubber sliding distance in the simulation is 114.3 mm, thus the change in temperature is negligible during the friction simulation in such a short time period of testing.
The surface-to-surface contact algorithm incorporated in the exponential decay friction model is used to simulate the surface interaction between the rubber slider and the pavement surface morphology. The resultant tangential force is computed for the rubber slider, where velocity and loading conditions are applied. The ratio of the tangential force to the normal load is defined as the simulated skid resistance [62], as shown in Figure 8. same compounding requirements as the standard tires used for common pavement skid-resistance testing devices, such as locked-wheel trailer and grip tester [41,[57][58][59]. However, rubber manufacturers usually do not publish material's mechanical property information. Many previous research studies have simulated the tire-pavement interaction by setting the rubber as a linear elastic material with a Poisson ratio around 0.5 [30,32,45,60]. It is also pointed out that the linear elastic material model still has the potential to demonstrate the rubber's mechanical behavior by [61]. Hence, in this research, the linear elastic material properties are adopted from the papers by Ong et al. and Zhang et al. [45,60], in which the tire rubber material property is as same as the DFT rubber slider's. The elastic modulus of the rubber slider is 100 MPa and the Poisson's ratio and density of the rubber are set as 0.45 and 1200 kg/m 3 , and rubber's viscoelastic parameter's is adopted from the papers by Zegard et al. [61], respectively. The applied load on the slider is 11.8 N, and the sliding speeds for simulation are 20, 40 and 60 km/h. Since all the data collection was completed with very similar temperature conditions, the properties of the rubber material are considered to be constant. The rubber sliding distance in the simulation is 114.3 mm, thus the change in temperature is negligible during the friction simulation in such a short time period of testing.
The surface-to-surface contact algorithm incorporated in the exponential decay friction model is used to simulate the surface interaction between the rubber slider and the pavement surface morphology. The resultant tangential force is computed for the rubber slider, where velocity and loading conditions are applied. The ratio of the tangential force to the normal load is defined as the simulated skid resistance [62], as shown in Figure 8. The water film thickness for the DFT testing is generally considered as 1 mm [41], thus the rough spots on the pavement surface with micro-rough surface characteristics are able to break-through the film of water present at the rubber-pavement interface and then form skid resistance [63]. In this research, the water's influence on pavement friction is considered by lowering the friction parameters in the solid-solid model as a thin lubricant in the rubber-pavement interface [43]. In the FEM simulation, the exponential decay friction is "lowered to represent the introduction of a lubricant between the bodies" [54]. Such assumption was originally proposed by Oden and Martins [34] and subsequently adopted in the ABAQUS ® software (Version 6.13, Dassault Systèmes, Vélizy-Villacoublay, Paris, France). Additionally, the derived exponential decay friction parameters are back-calculated from the in-situ collected DFT data, and thus the influence of water on the rubberasphalt (pavement) friction is included in the derived friction parameters.
Mesh study is performed to select the appropriate element size of the slider. Assuming the exponential decay friction parameters ( , , and α) to be 0.4, 0.35 and 0.6, the mesh study results in Figure 9 show the comparisons of the discretization errors and the computing times at the speeds of 20 km/h (DFT 20), 40 km/h (DFT 40), 60 km/h (DFT 60). It is found that the discretization errors keep less than 5% in the beginning with varying element sizes, and later significantly rises up to roughly 12% after element sizes reach 0.5 mm. The convergence could be achieved when the discretization error is less than 5% while the computing time remains low with various element sizes. As a result, considering both the numerical convergence and the computational efficiency, the rubber's element size is selected as 0.5 mm. The water film thickness for the DFT testing is generally considered as 1 mm [41], thus the rough spots on the pavement surface with micro-rough surface characteristics are able to break-through the film of water present at the rubber-pavement interface and then form skid resistance [63]. In this research, the water's influence on pavement friction is considered by lowering the friction parameters in the solid-solid model as a thin lubricant in the rubber-pavement interface [43]. In the FEM simulation, the exponential decay friction is "lowered to represent the introduction of a lubricant between the bodies" [54]. Such assumption was originally proposed by Oden and Martins [34] and subsequently adopted in the ABAQUS ® software (Version 6.13, Dassault Systèmes, Vélizy-Villacoublay, Paris, France). Additionally, the derived exponential decay friction parameters are back-calculated from the in-situ collected DFT data, and thus the influence of water on the rubber-asphalt (pavement) friction is included in the derived friction parameters.
Mesh study is performed to select the appropriate element size of the slider. Assuming the exponential decay friction parameters (µ s , µ k , and α) to be 0.4, 0.35 and 0.6, the mesh study results in Figure 9 show the comparisons of the discretization errors and the computing times at the speeds of 20 km/h (DFT 20), 40 km/h (DFT 40), 60 km/h (DFT 60). It is found that the discretization errors keep less than 5% in the beginning with varying element sizes, and later significantly rises up to roughly 12% after element sizes reach 0.5 mm. The convergence could be achieved when the discretization error is less than 5% while the computing time remains low with various element sizes. As a result,

Binary Search Back-calculation Approach
To obtain the exponential decay friction model's parameters µ s , µ k , and α, the binary search back-calculation methid is adopted in this study. This method consistently adjusts these parameters in the FEM simulation process so that the simulated skid resistance values are approximating the in-situ measurement value within acceptable accuracy. The DFT skid resistance measurements at 20 (DFT 20), 40 (DFT 40), and 60 km/h (DFT 60) are selected as the field validation data for this method. The limits of the 95% confidence interval are fulfilled when the percentage error from the back-calculation process is less than 10% [60].
Zhou et al. [33] recommended that the static friction coefficient µ s and the kinetic friction coefficient µ k to be 0.85 and 0.70 respectively, and the decay coefficient α ranging from 0 to 1 under dry conditions. Since the friction model in this paper incorporates the influence of water lubrication, µ s , µ k , and α should fall within the ranges of [0, 0.85], [0, 0.7], and [0, 1] [38].
Using the calculation of µ k as an example, the binary search back-calculation method is composed of the following steps [64]: (1) Calculate the midpoint value σ 0 from the initial range interval [0, 0.85]; (2) Input σ 0 into the FEM simulation process and obtain the simulated skid number SN σ 0 ; (3) Compare the SN σ 0 . with the in-situ skid number SN in−situ . If the percentage error (SN σ 0 − SN in−situ )/(SN in−situ ) is less than 10%, the convergence is satisfied and the iteration process stops. Otherwise the iteration process continues and proceeds to step 4; (4) Calculate the new midpoint σ 1 from the subinterval [0, σ 0 ]. If the difference from step 3 is positive or [σ 0 , 0.85] if the difference is negative, feed it into the FEM simulation process for simulated skid resistance SN σ 1 , and check whether the differences are within desired accuracy range. Repeat the process until the percentage error is less than 10%.
According to the exponential decay friction formula, µ k is the main factor affecting the interface friction at high speed, µ s dominantly affects the interface friction at low speed, the decay coefficient α represents the friction coefficient's sliding velocity dependency property. Therefore, it is computationally effective to firstly back-calculate µ k at the speed of 60 km/h. Subsequently, keeping µ k constant, back-calculate µ s at the speed of 20 km/h. Finally, keeping µ k , µ s constant, back-calculate α at the speed of 40 km/h, following the back-calculation process described above until the percentage error is less than 10%. For each data set, the back-calculation process is applied to obtain the appropriate exponential decay friction parameters via a loop computing process until the computed skid resistance best fit the in-situ DFT measurements. Table 2 provides one example of the back-calculated exponential decay friction model's parameters.

Validation of FEM Simulation Results
The DFT measurements at the speed of 20 km/h (DFT 20), 40 km/h (DFT 40), 60 km/h (DFT 60) are then used to back-calculate the exponential decay friction model parameters. As shown in Figure 10a, the R-squared values are 0.68 to 0.85 between the simulated and DFT measured skid resistance. To further validate the back-calculation method, the skid resistance at the speed of 30 km/h (DFT 30), 50 km/h (DFT 50), 70 km/h (DFT 70), which are not used for the back-calculation process, are simulated for verification. As shown in Figure 10b, the R-squared value varied from 0.55 to 0.68, indicating that the derived exponential decay friction model parameters are satisfactory for FEM simulation.

Validation of FEM Simulation Results
The DFT measurements at the speed of 20 km/h (DFT 20), 40 km/h (DFT 40), 60 km/h (DFT 60) are then used to back-calculate the exponential decay friction model parameters. As shown in Figure  10a, the R-squared values are 0.68 to 0.85 between the simulated and DFT measured skid resistance. To further validate the back-calculation method, the skid resistance at the speed of 30 km/h (DFT 30), 50 km/h (DFT 50), 70 km/h (DFT 70), which are not used for the back-calculation process, are simulated for verification. As shown in Figure 10b, the R-squared value varied from 0.55 to 0.68, indicating that the derived exponential decay friction model parameters are satisfactory for FEM simulation.

Model Development
Since the FEM process is time consuming and computationally extensive, it is desired to develop friction prediction models using in-situ 3D texture data sets and parameters. The detailed definitions of the 3D texture indicators at the macro-and micro-scales can be found in Table 3. These parameters have been widely used to characterize surface texture properties. The descriptive statistics of 3D areal pavement texture parameters are summarized in Table 4, which includes statistics such as average, maximum, minimum, standard deviation. Cross-correlation analysis in previous study [16] is conducted to reveal the correlation among the 3D macro-and micro-texture indicators within each category. It has been demonstrated that a high level of correlation exists within the macro-and micro-texture indicators. In order to enable accurate mapping of the texture indicators to the friction parameters, it is important to reduce the dimensionality of the identified 3D texture indicators and then develop a multivariate regression model for friction prediction. In this paper, principal component analysis (PCA) [65,66], a statistical procedure that uses an orthogonal transformation to convert a set of observations of possibly correlated variables into a set of values of linearly uncorrelated variables called principal components, is used for dimension reduction. The basic equation is defined as below: where, A is the PC vectors matrix; B is the original variable matrix; Q is weighting coefficients matrix to establish the linear relationship between A and B; m is the number of variables in original matrix; n is the number of experiments/observations. After the variables transformation, friction prediction models can be developed based on the significant PC vectors, as shown in Equation (3): where, a * is estimated coefficient for intercept; P * i is the PC vector; b * i is the corresponding coefficient for the PC vector.

Model Verification and Discussion
The derived parameters from the exponential decay friction model and the 3D texture indicators are used as the dependent variables and explanatory variables for the PCA regression process. 27 pairs of data sets are randomly selected to develop the regression model and 9 data sets are used to validate the model.
The PCA regression results are summarized in Table 5. The p-values for the PC vectors are all smaller than 0.01, implying that the generated principal components are significant to the friction parameters. The adjusted R square value of the regression model for µ k , µ s and α are 0.9456, 0.8276, 0.7215, respectively, demonstrating that the friction parameters can be well explained by the principal component vectors after eliminating the multicollinearity. Notes: Significant values: *, p < 0.5 × 10 −01 ; **, p < 0.1 × 10 −01 ; ***, p < 0.01 × 10 −01 . The more p-value with star marks, the higher significance exists between the principal component and exponential decay friction parameters.
To illustrate the influence of texture indicators on exponential decay friction parameters, the selected significant PC vectors are transformed back into the combinations of original texture indicators as shown in Table 6 so that the friction parameters can be directly predicted as shown in Equation (4): where, a is estimated coefficient for intercept; T is the texture indicators; b i is the corresponding coefficient for texture indicators. It can be concluded from Table 6 that texture indicators comprehensively influence the exponential decay friction parameters. Macro-texture and micro-texture commonly dominate the µ s , while micro-texture dominate the µ k as well as the macro-texture dominate the α. However, different texture indicator plays different roles on the friction parameters. For example, Sq dominate the µ s , µ k , and α, while Ssk significantly determine µ k at micro-level but fail to dominate the µ s and α at either micro-or macro-level. This phenomenon also implies that the individual texture indicator's influence on the friction parameter changes when switching from macro-level to micro-level.
The remaining data sets are used to compare the predicted values with actual DFT measurements, as shown in Figure 11. The differences between the predicted results and experimental data are mostly under 10%, which can be demonstrated by the largely overlapping of experimental and predicted values.
The predicted skid resistance has the similar trend with the in-situ skid resistance measured by the DFT at various testing speeds. It can be concluded that good agreements exist between the skid resistance predicted from the FE model in this research and the measurement results. This indicates that the proposed methodology can be used to predict skid resistance in terms of simulating the DFT with acceptable accuracy. The remaining data sets are used to compare the predicted values with actual DFT measurements, as shown in Figure 11. The differences between the predicted results and experimental data are mostly under 10%, which can be demonstrated by the largely overlapping of experimental and predicted values.   The predicted skid resistance has the similar trend with the in-situ skid resistance measured by the DFT at various testing speeds. It can be concluded that good agreements exist between the skid resistance predicted from the FE model in this research and the measurement results. This indicates that the proposed methodology can be used to predict skid resistance in terms of simulating the DFT with acceptable accuracy.

Conclusions
This study proposed a framework to quantify the relationship between texture characteristics

Conclusions
This study proposed a framework to quantify the relationship between texture characteristics and the interface friction coefficient and to prepare the friction inputs for the 3D based rubber pavement interaction simulation using testing data from the LTPP SPS-10 WMA testing site in Oklahoma. Comparing to the previous study [9,22], this research considers pavement texture as a significant factor for the skid resistance prediction at both macro-and micro-level. In particular, the following analyses are conducted, which could help researchers better investigate the rubber pavement interaction mechanism and aid road agencies making better pavement maintenance decisions: • A rubber areal pavement interaction FEM model is established to determine the rubber pavement interface friction by re-constructing 3D areal pavement model from high resolution surface texture data; • The binary search back-calculation method is used to derive the rubber pavement interface friction parameters so that the simulated skid resistance fits with the in-situ skid resistance data at a desired accuracy level; • PCA regression models are developed to correlate interface friction parameters and the 3D areal pavement texture characteristics, which can be used to prepare the inputs of friction parameters for FEM simulation.