A Spatial-statistical model to analyse historical rutting data

Rut depth is measured annually using the ViaPPS (Pavement Profile Scanner - laser scanner), developed by ViaTech in collaboration with the NPRA (Peraka and Biligiri 2020). Measurement positions are determined with high-accuracy GPS (CPOS/DPOS — satellite-based position service with centimetre/decimetre accuracy), ensuring precise location data. The rut depth is estimated by comparing the height of the ridges formed between the wheel paths, which involves identifying the two lowest points in the wheel path and measuring the distance to the top of the ridge. Additionally, the laser covers an approximate 4-meter footprint, ensuring that the vehicle remains within the wheel tracks and avoids crossing the centerline or shoulder. Thus, minor deviations (10-50 cm) do not affect the calculations. Annual rut depths along the road range from 0 to 41 mm, measured at 20-metre intervals. There are 3355 road segments measured for annual rut depth.

Rutting, denoted as $r_{i,t}$ (mm/year), is calculated as the difference in rut depth between years, as shown in Equation 1. We denote the rut depth variable at road segment $i=1,...,n$ in year $t$ as $d_{i,t}$, and $r_{i,t}$ as rutting (mm/year). Rutting for year $t=1,...,T$ at segment $i$ is defined as

When repavement occurs, rut depth decreases, sometimes resulting in large negative values of the rutting ($r_{i,t}$). Following Vedvik (2021), large negative rutting values are considered missing and set to NaN (not a number):

and the small negative values kept in the dataset are considered as measurement uncertainties and are accounted for in the random variability. Some missing rutting data may occur due to measurement uncertainty or data collection process. For available observations, the overall mean rutting and segment mean rutting (over all years) are calculated as $\bar{r}={1}/{n}\sum_{\forall t}\sum_{\forall i}r_{i,t}$ and $\bar{r}_{i}={1}/{T}\sum_{t=1}^{T}r_{i,t}$, respectively.

In the database, road width and the number of lanes are available. The first 2.5 km from Stjørdal has four (4) lanes, the remaining stretch has two lanes. Lane width is calculated by dividing the road by the number of lanes.

The Annual Average Daily Traffic (AADT) for each road segment between 2010 and 2020 varied from 568 to 18,500 vehicles per day (Figure 1, graph (a)). Fluctuations in traffic volume over short distances are influenced by various factors. Notably, there are considerably larger traffic volumes within Stjørdal to approximately 2.5 km towards Storlien. The European route E14 runs through the city center of Stjørdal, contributing to high traffic levels in the initial few kilometers.

The Norwegian Public Roads Administration (NPRA) recommends maintenance based on AADT values. High-traffic roads should be treated when rut depths exceed 20 mm, while low-traffic roads should be treated when rut depths exceed 25 mm. Adjusting AADT values based on road width is essential. Throughout the analyses, AADT values are expressed in ten-thousands.

Asphalt mixtures in cold-temperature regions like Norway are specifically designed to withstand harsh climatic conditions. These mixtures may include Polymer Modified Asphalt (PMA) to prevent cracking due to freeze-thaw cycles and low temperatures, and low-temperature performance grade binders to maintain flexibility at sub-zero temperatures. Dense-graded aggregates are used to minimize water infiltration, crucial for preventing freeze-thaw damage. The primary focus is on enhancing flexibility, low-temperature performance, and resistance to moisture, employing softer binders, dense-graded aggregates, and anti-stripping agents (Snilsberg et al. 2006, d’Angelo et al. 2008, Saba et al. 2009, Ødegård 2015).

Asphalt concrete (Ac), asphalt gravel concrete (Agc), and stone mastic asphalt (Sma) are the primary surface mixtures used in Norway. Sma, known for its superior wear resistance and deformation resistance due to the larger- and high-quality-stone in its composition, is typically employed on high-traffic roads. In contrast, Ac and Agc are commonly used on low-traffic roads, courtyards, sidewalks, and bicycle paths.

The model without spatial effects (non-spatial) has the linear predictor

The model assumes that the level of expected rutting (mm/year) is influenced by an interaction effect of the annual average daily traffic (AADT) and some explanatory variables, including the level of rut depth from the previous year ($z_{\mathrm{d-1}}$), the lane width ($z_{\mathrm{w}}$) and asphalt type. This assumption stems from our prior understanding: when there’s no traffic, there’s no rutting (or rut depth). This insight is grounded in research, such as the work of Hao et al. (2009), which identified increased loading time (traffic loads) as the primary factor accelerating rutting in asphalt pavements.

Let $z_{\mathrm{AADT}(s,t)}$ denote AADT for location $s$ in year $t$, and $I_{\mathrm{Ac}(s,t)}$ be an indicator variable that has a value of 1 if location $s$ in year $t$ has asphalt type Ac and $0$ otherwise. Then, $z_{\mathrm{Ac}(s,t)}=I_{\mathrm{Ac}(s,t)}\cdot z_{\mathrm{AADT}(s,t)}$ and, hence, equals AADT if the asphalt type is Ac (asphalt concrete) and $0$ otherwise. Similarly, $z_{\mathrm{Agc}(s,t)}$ and $z_{\mathrm{Sma}(s,t)}$ equal the AADT for location $s$ in year $t$ if the asphalt type is Agc (asphalt gravel concrete) and Sma (stone mastic asphalt), respectively, and $0$ otherwise. Further, $z_{\mathrm{d-1}(s,t)}=d(s,t)\cdot z_{\mathrm{AADT}(s,t)}$ and $z_{\mathrm{w(s,t)}}=\mathrm{w}{(s,t)}\cdot z_{\mathrm{AADT}(s,t)}$.

The random yearly effect is assumed to be distributed as $\displaystyle\boldsymbol{\gamma}=(\gamma_{1},\gamma_{2},...,\gamma_{T}){\sim}\mathrm{Normal}(\boldsymbol{0}_{T},\ \tau_{\gamma}^{-1})$.

Our model (Equation (3)) does not include several variables known to influence rutting, such as weather, bearing capacity, and proportion of heavy traffic load. These variables are spatially dependent, with segments closer in space likely to be more similar than those further apart. To account for this spatial dependency, we incorporate a Gaussian process (GP) as a spatial random effect along the road. GPs are powerful tools for modeling spatial dependencies, capturing variability due to both known factors (e.g., road characteristics, traffic patterns, and other explanatory variables) and unknown factors not explicitly accounted for in the model. The unknown factors could include unmeasured road characteristics, environmental variables like temperature, or other spatially dependent factors influencing rutting. Our model aims to capture both consistent spatial effects that remain stable over time and spatial effects that vary between years.

Consistent spatial effects, such as road characteristics (e.g., bearing capacity, side slope height, surface and base curvature, and proportion of heavy traffic), are likely to influence rutting consistently over time. Additionally, there are spatial effects that vary between years, such as heavy rain events, temperature variations, or fluctuations in heavy traffic due to construction work. To account for these variations, we include both an annually varying spatial process, $\xi_{t}(\boldsymbol{s})$, for each year $t=1,...,T$, and a common spatial process, $\omega(\boldsymbol{s})$, referred to as the spatial rutting component.

The annually varying spatial process, $\xi_{t}(\boldsymbol{s})$, captures the spatial variability in rutting that changes from year to year, accounting for spatial variations specific to each year, such as the impact of specific weather events or construction activities.

The common spatial process, $\omega(\boldsymbol{s})$, which we refer to as the spatial rutting component, captures the spatial variability in rutting that remains consistent across all years. This includes spatial effects not explicitly accounted for by other variables in the model, such as road characteristics and other spatially dependent factors. This contribution can be either positive or negative. The linear predictor for the model including spatial effects is described as

where the remaining model parameters are defined as that in Equation 3.

The common spatial process $\omega(\boldsymbol{s})$ in Equations (4) is the time constant variable. It is a Gaussian random process $\left\{\omega(\boldsymbol{s}):\boldsymbol{s}\in\mathcal{D}\subseteq\mathbb{R}^{d}\right\}$ such that for any $n\geq 1$ and for each set of spatial locations $\left(s_{1},...,s_{n}\right)$ satisfies

where $\boldsymbol{\mu}^{\omega}=\left\{\mu^{\omega}(s_{1}),...,\mu^{\omega}(s_{n})\right\}$ is the mean vector and $\Sigma^{\omega}_{ij}=\mathrm{Cov}\left\{\omega(s_{i}),\omega(s_{j})\right\}=C^{\omega}\left\{\omega(s_{i}),\omega(s_{j})\right\}$ are the elements of the covariance matrix defined by the Matérn stationary isotropic covariance function. Simply, this covariance is expressed as $C^{\omega}\left(d;\sigma^{2}_{\omega},\rho_{\omega}\right)$, where $d=||s_{i}-s_{j}||$ is the distance between two road segments; $\sigma^{2}_{\omega(s)}$ is the marginal variance of the spatial rutting component and $\rho_{\omega}$ is the spatial range for the spatial rutting component. These are the parameters to be estimated.

Similarly for the annually varying spatial field $\xi_{t}(\boldsymbol{s})$, the parameters to be estimated are $\sigma^{2}_{\xi_{t}}$ and $\rho_{t}$. Detail description of the Matérn covariance structure is found in (e.g., Blangiardo and Cameletti 2015, Krainski et al. 2018) and the Supplementary Material LABEL:supplementary:modelformula_SUPPLEMENTARY.

We consider simpler non-spatial and spatial models for comparison. These models include some of the known explanatory variables and are presented in Table 2. The most complex spatial model (in Equation 4) we refer to as Model 1, and it’s non-spatial version, Equation (3), as Model 4. These models include all three explanatory variables.

Asphalt types, particularly those with weaker binder concentrations like asphalt gravel concrete ($\hat{\beta}_{\mathrm{Agc}}$) and asphalt concrete ($\hat{\beta}_{\mathrm{Ac}}$), experience more rutting (Table 3). Qian et al. (2019) and Mascarenhas et al. (2020) shed light on how temperature and heavy traffic ($>$ 3.5 tonnes) impact rutting resistance on asphalt types such as Agc and Ac. Notably, asphalt gravel concrete is predominantly situated approximately 30 kilometers from the Swedish border, where a sub-arctic climate, underdeveloped networks, and heavy traffic converge to adversely affect pavement performance Peel et al. (2007), Ketzler et al. (2021), Jenelius (2010), Hovi et al. (2018), Pinchasik et al. (2020). There is also greater uncertainty in asphalt concrete ($\hat{\beta}_{\mathrm{Ac}}$) and asphalt gravel concrete ($\hat{\beta}_{\mathrm{Agc}}$) compared with stone mastic asphalt ($\hat{\beta}_{\mathrm{Sma}}$), although considering spatial dependencies helps to explain more of the variability. Non-spatial models appear to explain variability between years rather than by asphalt type, suggesting a nuanced relationship between these factors in understanding rutting susceptibility.

The rut depth from the previous year ($\hat{\beta}_{\mathrm{d-1}}$) becomes increasingly important when spatial variability is not considered (Model 4, Table 3). As a genuine spatial-temporal variable, including the previous year’s rut depth in the model modifies the spatial effect $\hat{\sigma}_{\omega}(s)$, specifically reducing it, which in turn enhances the precision of the estimates by decreasing uncertainty, $\hat{\sigma}_{\varepsilon}$. This improvement is evident in Models 1 and 3 (see Table 4 in Section A). Incorporating additional relevant spatial-temporal variables could further refine the results, providing a better explanation for the observed uncertainty.

The estimated effect of lane width ($\hat{\beta}_{\mathrm{w}}$) is minimal and positive, suggesting negligible impact on expected rutting.

The observed mean rutting (Figure 2 (a) in Section 2.3) stands at $1.62$ mm/year. In Figure 4, graph (a), 361 road segments have an expected rutting exceeding this mean, with 49 segments showing estimated rutting exceeding $3$ mm/year (graphs (c) and (d), respectively). These 49 locations signify potential hot spots for accelerated rutting, necessitating urgent maintenance, as their expected lifetime could be as short as 8 years, depending on traffic. Additionally, the observed mean rutting across segments ranges from -0.14 to 8.33 mm/year (Figure 2, graph (a)), while the maximum expected rutting from the model is 4.67 mm/year (Figure 4, graph (a)), approximately half the maximum observed mean rutting.

In Figure 4, graph (b), we plot the expected rutting $\hat{\eta}$, lane width, and AADT (per ten thousand) for a short road section spanning road stretch 42000 to 54000 metres (or $42-54$ km from Stjørdal), characterized by weaker asphalt type Agc. This section exhibits spikes of accelerated rutting, particularly evident with increasing traffic. Possible factors contributing to this variability include changes in asphalt type, traffic composition, suboptimal road design, or inadequate pavement maintenance.

To further investigate high rutting values in a road section, we analyze the impact of explanatory variables, yearly variation, and spatial rutting component ($\hat{\eta}$). We explore spatial Model 1 using different values for variables to estimate rutting. First, we use actual measurements like asphalt type, lane width, and rut depth ($\hat{\eta}_{\mathrm{True}}$). Second, we consider alternative asphalt types - asphalt concrete (Ac) or stone mastic asphalt (Sma) ($\hat{\eta}_{\mathrm{Ac}}$ or $\hat{\eta}_{\mathrm{Sma}}$). Third, we examine extreme values such as the more wear- and deformation-resistant pavement (Sma), minimal rut depth (0.01 mm), or wider lanes (6 meters) ($\hat{\eta}_{\mathrm{Best}}$). Conversely, accelerated rutting scenarios involve deeper rut depth (25 mm), narrower lanes (2.75 meters), and asphalt gravel concrete (Agc) ($\hat{\eta}_{\mathrm{Worse}}$) from the E14 dataset. These scenarios are depicted in Figure 5 (a).

With worse values ( ), rutting is projected to be at least 1.5 mm/year ($1.76-8.20$), indicating accelerated deterioration and the need for maintenance within 3 to 14 years. Locations with unexpectedly high rutting, exceeding 4 mm/year, are observed, possibly due to increased AADT compared to other areas (Figure 4, graph b).

Using the best covariate values ( ), rutting is projected to be at most 4.30 mm/year, with maintenance required approximately 6 years after construction or maintenance. Additionally, using the asphalt binder Sma yields similar rutting estimates but generally lower than those using the best values, indicating asphalt type as a primary factor. Asphalt gravel concrete (Agc,  ) and asphalt concrete (Ac,  ) show similar rutting estimates due to similar binder concentrations.

We examine the impact of asphalt type on rutting (Figure 5, graph (b)). Changing from asphalt type Agc (asphalt gravel concrete) to the more wear- and deformation-resistant Sma (stone mastic asphalt) while keeping all other variables constant, we sequentially add additional explanatory variables to the model, starting with asphalt type. Initially, asphalt type (Agc or Sma) is included with the random yearly effect $\hat{\gamma}$ and the spatial rutting component $\hat{\omega}(s)$ as follows: $\hat{\tilde{\eta}}_{\mathrm{Agc}}=\hat{\beta}_{\mathrm{Agc}}z_{\mathrm{Agc(s,t)}}+\hat{\gamma}+\hat{\omega}(s)$ or $\hat{\tilde{\eta}}_{\mathrm{Sma}}=\hat{\beta}_{\mathrm{Sma}}z_{\mathrm{Agc(s,t)}}+\hat{\gamma}+\hat{\omega}(s)$. Subsequently, the previous year’s rut depth is incorporated, represented as $\hat{\tilde{\eta}}_{\mathrm{d-1}}=\hat{\beta}_{\mathrm{Agc}}z_{\mathrm{Agc(s,t)}}+\hat{\beta}_{\mathrm{d-1}}z_{\mathrm{d-1(s,t)}}+\hat{\gamma}+\hat{\omega}(s)$ or $\hat{\tilde{\eta}}_{\mathrm{d-1}}=\hat{\beta}_{\mathrm{Sma}}z_{\mathrm{Sma(s,t)}}+\hat{\beta}_{\mathrm{d-1}}z_{\mathrm{d-1(s,t)}}+\hat{\gamma}+\hat{\omega}(s)$. Similarly, the model with added lane width ${\hat{\beta}_{w}z_{\mathrm{w(s,t)}}}$ is denoted as $\hat{\tilde{\eta}}_{\mathrm{w}}$. Results from using stone mastic asphalt (Sma) are depicted by solid lines, and asphalt gravel concrete (Agc) results are shown with dotted lines.

The findings affirm that asphalt type significantly influences rutting. Implementing the binder type Sma instead of the actual asphalt type Agc would decrease expected rutting, reaching a maximum of $4.35$ mm/year (compared to 4.68 mm/year), approximately reducing rutting by 0.33 mm/year annually. Over a projected lifetime of 6 years, this equates to an additional year of serviceability.

We examine bearing capacity, surface curvature index (SCI, micrometre), and base curvature index (BCI) for road segments 42000-54000 using 2021 data, indicative of construction quality and soil pressure sustainability over the past decade. Bearing capacity is a measure of the deflection caused by a rolling load axle on pavement. It is measured using the Rapid Pavement Tester (Raptor, https://www.ramboll.com/rst/rst-products), which continuously collects data on the pavement structure. A value exceeding 16 metric tons (       ) suggests a robust foundation design (Figure 5, graph (c)). The SCI to BCI ratio (Figure 5, graphs (d) and (e)) reveals deflection location within the pavement layers. High rutting may result from factors like side slope height (data unavailable), inadequate drainage, or sub-base layer deformation during thawing periods. Bearing capacity data, sampled and temperature-corrected during summer, may inflate measurements compared to thawing periods. Varying construction material quality along the road stretch (30000-67100 metres) warrants investigation into surface and subsoil layer deflections. An SCI/BCI ratio below 12 (       ) indicates robust construction. Inconsistent ratios suggest layer deformation, with data indicating uneven construction on E14, featuring deeper base layer deflections.

In conclusion, asphalt type significantly influences rutting, recommending more wear resistant and deformation resistant binder types for extended lifespan.