ROAD TEMPERATURE MODELLING WITHOUT IN-SITU SENSORS

Modelling of the pavement temperature facilitates winter road maintenance. It is used for predicting the glaze formation and for scheduling the spraying of the de-icing brine. The road weather is commonly forecasted by solving the energy balance equations. It requires setting the initial vertical profile of the pavement temperature, which is often obtained from the Road Weather Information Stations. The paper proposes the use of average air temperature from seven preceding days as a pseudo-observation of the subsurface temperature. Next, the road weather model is run with a few days offset. It first uses the recent, historical weather data and then the available forecasts. This approach exploits the fact that the energy balance models tend to “forget” their initial conditions and converge to the baseline solution. The experimental verification was conducted using the Model of the Environment and Temperature of Roads and the data from a road weather station in Warsaw over a period of two years. The additional forecast error introduced by the proposed pseudo-observational initialization averages 1.2 °C in the first prediction hour and then decreases in time. The paper also discusses the use of Digital Surface Models to take into account the shading effects, which are an essential source of forecast errors in urban areas. Limiting the use of in-situ sensors opens a perspective for an economical, largescale implementation of road meteorological models.


Introduction
Winter weather conditions significantly influence the road safety. The most dangerous phenomena, such as the formation of glaze or ice, are counteracted by ploughing and spraying actions. Their scheduling sometimes is supported by specialized weather forecasts for roads. Although the pavement temperature is of central interest, the traffic safety depends also on other factors such as the presence of water on the road surface.
The forecasting of the pavement temperature with physical models requires data from the Road Weather Information Station. They are used for the initialization of the vertical temperature profile and adjusting the weather forecasts to the local conditions. This paper investigates an alternative approach, in which only the widely available weather forecasts and recent historical data are used. This study provides an analysis of the temperature convergence speeds and the observed biases. It also suggests the seven-day average air temperature as a reasonable estimate for the initial subsurface pavement temperature. The proposed methods were demonstrated in urban areas (cities Bordeaux and Warsaw), where the shades from buildings and trees were a significant source of prediction imprecision. The paper shows how they are taken into account by appropriate processing of the Digital Surface Models.

Road weather forecasting system METRo
This study is based on the Model of the Environment and Temperature of Roads (METRo), version 3.3.0. Environment Canada distributes is as an open source software under the GNU General Public License, version 2. METRo has been developed from 1999 and is being introduced in other countries cf. (Karanko et al. 2012;Kršmanc et al. 2012;Linden, Petty 2008;Sokol et al. 2014). Crevier and Delage (2001) documented the physical and implementation principles in an article. Model of the Environment and Temperature of Roads computes the energy balance of the road surface taking into account energy fluxes schematically depicted in Fig. 1. More specifically, the atmospheric side is modelled using the following Eq (1): where R -residue, W·m −2 ; α -albedo, 1; S -incoming solar flux, W·m −2 ; ε -emissivity, 1; I -incoming infrared radiation, W·m −2 ; σ -Stefan-Boltzmann constant, W·m −2 ·K −4 ; T s -pavement surface temperature, K; H -sensible turbulent heat, W·m −2 ; L a -vaporization heat, J·kg −1 ; E -water vapour flux, kg·m −2 s −1 ; L f -heat of fusion, J·kg −1 ; P -water precipitation rate, kg·m −2 s −1 ; A -anthropogenic flux, W·m −2 . Residue R is computed by summing the net solar flux with the net infrared radiation, obtained as the difference between the absorbed incoming flux εI and the emitted flux −εσ 4 s T . Next, the sensible turbulent heat flux H and the latent heat flux are subtracted. The flux ± f L P corresponds to the phase changes of water precipitating with rate P, where the sign ± denotes freezing or thawing. Anthropogenic flux A originates, e.g., from passing cars (Crevier, Delage 2001).
METRo also contains a module describing heat conduction within the road layers using one-dimensional heat diffusion Eq (2): where C -the volumetric heat capacity, J·m −3 ·K −1 ; zdepth, m; T -temperature, K; t -time, s; G -ground heat flux, W·m −2 . Equation (3) gives the ground heat flux G: where k -the heat conductivity, W·m −1 ·K −1 .
Since constants C and k are specific to different road materials, they change at depths z corresponding to different road layers (Crevier, Delage 2001). The upper boundary condition is the rest downward flux R. However if freezing or thawing is taking place, it is changed to a requirement that surface temperature equals 0°C (273.15 K). The bottom boundary condition is zero energy flux for normal roads or air temperature T a , K for bridges and overpasses.
The METRo model also contains a module for water accumulation. Two compartments depict the amounts of the accumulated water in liquid and solid state. This approach is based on the results of Saas (1992).

Convergence of road temperature forecasts
Running a model requires specifying the initial conditions for the differential equations describing the accumulation and conduction of energy. METRo uses a roadside meteorological station to obtain the initial temperature of pavement surface and subsurface at a depth of 40 cm. These data are used for the initialization of the vertical temperature profile and for coupling the observations with forecasts. If more than four hours are missing, an analytical solution to the heat conduction equation is used based on the diurnal cycle of fixed amplitude. Crevier, Delage (2001) note: "Results using this analytic road temperature profile were surprisingly close to those using the full initialization method. The difference in road temperature after a 24 h forecast was never greater than 0.18 °C in tests that compared forecasts with full initialization to ones for which no initialization data were supplied".
The data from a road meteorological station increase the accuracy of the forecasts and allow for adjustments to the local conditions. However, running road weather forecast system without in-situ data would bring clear benefits: − no costs of building and maintaining additional meteorological infrastructure; − the possibility to forecast road weather for large networks. The current and historical data are readily available from many meteorological stations, in particular, those located at the airports. Similarly, numerical weather forecasts are provided by professional services. They serve as a substitute of the road weather stations. The only missing parameters are the surface and subsurface temperatures. However, appropriate processing of the historical data is a way to approximate them.
In Fig. 2, nine surface temperature forecasts are shown for the same period in winter (January 12 th -18 th ). They differ in the choice of the initial surface and subsurface temperatures, which ranged from −20 °C up to +20 °C by 5 °C. The maximal difference among the forecasts declined in time, form the initial +40 °C to +13.5 °C, +8.6 °C, and +3.1 °C after one, two, and three days. On the fourth day, the difference increased to +5.5 °C, as the weather was The example in Fig. 2 suggests that the model "forgets" its initial conditions. For instance, the overestimation of the surface temperature would increase the emission flux −εσT s 4 strengthening the cooling flux in the forecasts. Consequently, the predictions converge to the oscillatory equilibrium. Formally, this statement is derivable using the Lyapunov Second Method for stability applied to the system of differential Eqs (2)-(3). However, there is also a clear physical interpretation. Imagine that the initial temperature was inflated over its (dynamic) equilibrium through an anthropological action such as a leakage of a hot fluid from a heat network. The road is expected to cool down, what is observed in Fig. 2. Analogous reasoning applies to the case of underestimated temperature. Consequently, the in-situ measurements are not indispensable. Prognosis is possible for any initial temperature, if the forecast is run for a few days before the current date. Bouilloud et al. (2009) successfully applied a similar simulation-based approach.
In the default setting, METRo uses only a limited amount of historical data, which overlaps with the last available weather forecast. Typically, this accounts for a few hours of overlap in the coupling stage (Crevier, Delage 2001;Linden, Drobot 2010). It allows indirect taking into account some of the local features of a given location (Karsisto et al. 2016). In the setting proposed in this paper, historical weather data from a few last days is used to estimate the plausible initial conditions. In Fig. 3 the data requirements are compared on a time scale.
A simulation experiment was run to quantify the rate of convergence of the pavement temperatures. It was based on historical weather data from winters 2013/14 and 2014/15 from the Warsaw Chopin Airport (ICAO code: EPWA). For each week, the model was run twice with the initial temperatures of surface and subsurface set to 5 °C below and 5 °C above the current air temperature. The convergence is illustrated in Fig. 4. The surface temperatures are shown with solid orange lines. A thick, dark orange, solid line indicates their mean. The yellow belt shows the band covering two standard deviations of the differences (mean ± standard deviation) indicating the level of variability of the data. During the first 24 hours, the difference among the surface temperatures decreases from +10 °C to +3 °C. Next, it gradually steps down showing oscillatory behaviour with a daily frequency. The oscillations result from a composition of a slowly changing heat conduction flux and the external fluxes, which strongly depend on the time of the day. Consequently, during days the temperatures converge, while in nights they slightly diverge. After 7 days, the average difference is less than 1 °C.
The subsurface temperatures at a depth of 40 cm are denoted in Fig. 4 by dashed lines. Their convergence process is much slower. After 7 days, the temperature difference is still above +5 °C. In the METRo model, the heat conduction is the only possibility of underground transferring of energy. The fact that the difference in surface temperatures levels off at non-zero value is related to the prolonged cooling (or heating) times by the underground road layers. The underground temperature curves are also much smoother, what indicates that the cooling process averages out the surface variations.
The sources of historical data are the in-situ measurements, a nearby road weather station, or a numerical weather model. Spatial extrapolation decreases the accuracy of results but allows for the easy and inexpensive application of the system to large road networks. Thermal mapping (Shao et al. 1997;Ząbczyk 2006) is a way to improve the forecast accuracy by taking into account the local variations.

Empirical verification
To thoroughly verify the pseudo-observational initialization, the results were compared against measurements obtained from a road weather station located at Aleja Krakowska, one of the main streets in Warsaw. Its exact location is shown in a screenshot in Fig. 5 from a system providing road data. The atmospheric weather observations originated from the nearby airport. Air temperatures measured in both stations were used to double-check the data and align timestamps. In a few cases, differences of up to 20 minutes were observed.
The METRo model is initialized with the surface and subsurface pavement temperature. The latter indicates the amount of thermal energy stored in the road structure. If it is incorrectly set then an additional, artificial heat flux is introduced into the model. Setting the subsurface temperature possibly close to its actual value minimizes this bias. The average air temperature from a preceding period is a substitute of subsurface sensors, which are often unavailable. Taking multiples of 24 hours ensures independence from the daily weather periodicity. For the analysed station the optimal number of days is n = 7. It was obtained by choosing the maximal Pearson correlation coefficient between the computed average and the subsurface temperatures, as shown in Table 1. In the analysed experimental setting, the subsurface sensor was placed at a depth of 30 cm, what is 10 cm shallower than the default setting in the METRo model.
The average air temperature from the preceding week was used as a pseudo-observation for the pavement surface and subsurface temperatures. A comparison of the measured surface temperature with its forecasts is shown in Fig. 6. The influence of the initialization method quickly decreases with time. After several hours, the forecasts practically overlap. The remaining discrepancy between the predictions and observations are due to other reasons.
Imperfect initialization is one of the sources of the prediction error. A simulation experiment was run for the winter periods of 2013-2014 using the data from the discussed road weather station to quantify its influence. Figure 7 compares the distribution of errors obtained for different initializations. The thick lines show the mean errors, whereas the bands cover two standard deviations of the errors obtained for each time instance. The initialization with observations from the road weather station gives slightly better results, but the overall accuracy is comparable, as other sources of imprecision are dominant.
The scheduling of the spraying and ploughing actions is usually based on short-term forecasts. The first few hours  are the most important for winter maintenance. Therefore, it is interesting to observe the influence of the approximate initialization in time Eq (4). Let e obs. (t) denote the difference between the surface temperature at time t and its METRo forecast obtained after initialization using observed data. Similarly, let e pseudo (t) denote the prediction error obtained with pseudo-observational initialization. The additional error e a (t) introduced by the approximate initialization is a difference between the absolute values of the errors: The additional error was estimated using the data from winter periods of 2013-2014 and is given in Table 2. Initially, it averages 1.2 °C and then steadily decreases in time. This result shows the expected accuracy loss due to using only the atmospheric weather data rather than the road information stations.
A daily periodic pattern is apparently visible in Fig. 7. The most substantial errors are observed during the daytime when direct insolation is the primary energy flux, and the temperatures change fast. The forecasts tend to be slightly out of phase (Linden, Drobot 2010), what further increases the daytime error. The underestimated anthropogenic flux is another source of variability in the presented results, as the analysed roadside road weather station is located at an arterial road next to a pedestrian crossing. Nevertheless, the observed error is quite repetitive, and subtracting its trend would improve the prediction accuracy. Introducing such corrections would change the nature of the model from purely physical to a mix of physical and statistical approaches. Such a change opens new perspectives, such as the estimation of confidence intervals for the forecasts (Sokol et al. 2017). Noteworthy, purely statistical methods provide comparably accurate results (Kršmanc et al. 2013). However, they require data from at least one harsh winter to fit the model.

Shading effects and water accumulation
Sky closure by roadside trees, buildings, or mountains affects the local pavement temperature through reducing the direct insolation. Neglecting this phenomenon can lead to significant prediction bias, especially in winter, when the sun is low in the sky. Kršmanc et al. (2014) upgraded the METRo model by introducing the shading effects. The algorithm is based on the outline of the visible horizon. For each azimuth an elevation is provided, above which the sky is clear. Consequently, the direct solar flux is accounted for only when the sun is above the visible horizon. Geographical location and time zone are used for the computation of the apparent motion of the sun across the sky.
The shading effects are particularly significant in urban areas, which were the test sites analysed in this study. To introduce them the elevations were computed for each azimuth at locations spread along the axes of the analysed roads. For this purpose, the Digital Surface Model (DSM) was used. It contains information about the heights of buildings, roads, trees and other structures. Fig. 8 compares an ortophotomap of a centre of Bordeaux with a corresponding DSM. Coarse scaled DSMs are available from satellite imagery. Sentinel 1 operated by the European Space Agency within Copernicus Programme provides data free of charge. More accurate images are provided by other missions such as GeoEye-1 or WorldView-2 (Aguilar et al. 2014). Fine-scaled surface models (with a pixel size of 1 m 2 or smaller) are available from airborne measurements cf. (Sirmacek et al. 2012). High-resolution DSM data is already available in many countries (Hu et al.  2016) and allows for detailed analysis of the shading and reflection effects. For a given azimuth ϕ i , the elevation was found as the maximal angle above that the sky is visible in a 1° cone within a range from R min = 5 m to R max = 50 m. In the polar coordinate system with the origin at the location of the observer, the analysed cone reads: where radius r is the distance from the observer. The set in Eq (5) is schematically illustrated in Fig. 9. The lower bound was introduced to reduce the influence of noise among individual pixels. The upper bound was set to 50 meters. The difference between the road weather forecasts computed with and without shades is exemplified in Fig. 10. The pavement temperature is generally cooler and more rugged in the shaded case. The surface temperature is also more likely to fall below the dew point leading to condensation of water vapour and changing the overall road condition from dry to wet.
Water has high thermal capacity, therefore its accumulation influences the pavement temperature. It also directly increases traffic hazards due to the possibility of aquaplaning or ice formation. Localization of places, where water accumulates on the road is often based on transversal cross-sections of the pavement, compare e.g. German standard ZTV ZEB-StB 2006 Zusätzliche Technische Vertragsbedingungen und Richtlinien zur Zustandserfassung und -bewertung von Straßen, in German. This method is only an approximation, therefore the resulting parameter is called Fiktive Wassertiefe (ZTV ZEB-StB 2006), i.e., fictional water depth. A two-dimensional height map of the pavement surface was assembled to localize water accumulation areas more accurately. It was based on longitudinal and transversal evenness measured by laser detectors mounted on the measurement vehicles.
Gyroscopes mounted in the vehicles measuring the pavement evenness also provide the inclination along the road. These data require special care since even slight calibration error accumulates along the route. It was corrected with the use of a DSM a reference source for the overall road height pattern. An alternative method is to correct the gyroscope bias by making a roundtrip (e.g., measuring road in both directions) and equate the accumulated height change to zero. The pavement surface model was then processed with the filling-transform algorithm. Fig. 11 shows that the area of water accumulation obtained in this way is much smaller than the "fictional" area calculated from transversal crosssections. It is also smoother and to a greater extend resembles paddles observed in reality. This method has a potential to improve the choice of the maintenance sections, especially in hilly areas where water flows along the road leading down. Nevertheless, issues such as the collection of water by kerbs in case of ineffective (e.g., clogged) drainage systems or water removal by passing cars remain for further studies.

Conclusions
1. The prediction of the pavement temperature is possible without roadside meteorological stations. Instead, averaged historical weather data from a few last days serve as the pseudo-observations of initial temperatures. This approach exploits the property that the energy balance models "forget" their initial condition.
2. Modelling the pavement temperature without insitu sensors is a very economical way of road weather forecasting for large networks, however, this comes at the cost of accuracy loss.
3. Average air temperature from the preceding 7 days is an adequate pseudo-observation for the surface and subsurface temperatures of the pavement.
4. The use of pseudo-observations initially increases the forecast error by an average of 1.2 °C with standard deviation 1.7 °C. In the following hours, this additional error decreases.   . 11. Comparison of water accumulation places obtained from cross-sections (above) and 2D pavement surface models (below) 5. The forecasts in urban or hilly areas are significantly improved by taking into account shading effects computed using the Digital Surface Models.