Exploring the Highway Icing Risk: Considering the Dynamic Dependence of Icing-Inducing Factors

： In terms of the dynamic dependence between icing-inducing factors, this study is to explore the risk 8 distribution of highways when icing events occur in the study area. A joint distribution considering the dynamic 9 correlation of inducing factors was first constructed employing the Copula theory, which then yielded the 10 possibility of icing events. Meanwhile, hazard zones and intensities of icing were proposed under different 11 exceeding probabilities. After finishing the vulnerability analysis of highways, the risk matrix was used to conduct 12 the icing risk for the highway, which was then applied to the construction of the risk zoning map. The results 13 showed that there was an upper-tail dependence between extreme precipitation and temperature in the study area in 14 winter, which could be well captured by the Gumbel Copula function. Indeed, the constructed joint distribution can 15 express the possibility of icing under different intensities of precipitation and temperature. Besides, the highway 16 with the tallest vulnerability in the study area was the Hegang-Yichun line. The case application showed that during 17 March 2020, the traffic lines with a high icing risk were distributed around Fujin, Jiamusi, Hegang, and Qitaihe 18 cities, and the Hegang section of the Hegang-Yichun line was at the highest icing risk. The low-risk lines were 19 concentrated in the western part of the study area. This study is of great significance for the prevention and control 20 of ice-snow disasters on the highway in cold regions.


Introduction
Precipitation causes water or snow on the road surface, which provides material for road icing. When the temperature 71 meets the icing conditions, the road surface with less precipitation will only maintain thin ice. As the precipitation 72 increases, the road surface develops in the direction of ice coverage (Figs. 1e). Given the data availability, the average 73 monthly precipitation (AMP) was used as an evaluation indicator for the material basis when road icing occurs. 74 Pavement temperature is the controlling factor of road icing. Together with precipitation, it determines the 75 development of road states, such as dryness, dampness, stagnant water and snow, slush, partial icing, and final ice 76 cover (Crevier and Delage, 2001;Berrocal et al., 2010). The temperature of the road surface is determined by the 77 temperature of the atmosphere. Here, the monthly average atmospheric temperature (AAT) was used as another 78 indicator (Nantasai and Nassiri, 2017). 79 As the traffic lines are exposed to the natural environment, the operation of highways will be affected when an icing 80 event occurs, such as vehicle speed reduction, road congestion, etc (Berdica, 2002;Somayeh et al., 2015). In this 81 study, highways in different locations of the study area were used as disaster-bearing objects (Fig. 1c). For 82 disaster-bearing bodies with different vulnerability characteristics, different risks will arise under the interference of 83 icing with the same intensity (Birkmann, 2015). Compared with the less vulnerable disaster-bearing body, the more 84 vulnerable hazard-bearing body is at a higher risk when involving the same hazard. Therefore, the vulnerability 85 analysis of highways is an indispensable step in the study of road icing risk. 86 Based on the disaster risk theory, the risk caused by the icing event to the highway in the study area can be divided 89 into two parts, including the hazard and the vulnerability analysis (Beven et al., 2015). In this study, AMP and AAT 90 were regarded as the two main inducing factors for road icing. The joint distribution function of the inducing factors 91 was established via the Copula theory, which is then used for the hazard analysis (Schölzel and P, 2008;Feng et al., 92 2017 The joint distribution function is a probability function, which describes the probability that the dependent variable 95 occurs under the combined action of multiple independent variables (Nelsen, 1997;Hotta, 2006;Hu, 2006

Methods
When using the MSED to determine the optimal Copula function, the Kendall rank correlation coefficient and the distance value was determined as the optimal Copula function to characterize the dependence between variables. 135 (2) 136 where 2 ( , ) n k d C C  represents the square Euclidean distance between n C  and k C , n C  depicts the empirical Copula, n is 137 the number of samples, k C depicts the candidate Copula function, u and v are the edge distribution function of 138 variables x and y, respectively, x and y represent the AMP and AAT. 139 Copula theory shows that any n-ary joint distribution function can be decomposed into n corresponding unary 141 distribution and a connection function describing the dependence between variables. In this study, the joint 142 Since the hazard is the product of abnormal changes in the natural environment, the hazard analysis and hazard 148 zoning of icing events must be related to the spatial differentiation of the natural environment (Fotheringham et al.,149 2017; Kumar et al., 2017). In this study, the hazard analysis of the regional icing events was conducted based on the 150 comprehensive zoning principle and the factor-dominant zoning principle of the territorial differentiation law. 151 In the principle of comprehensive zoning, different hazard zones were carried out based on the exceeding probability 152 of icing events (Zscheischler and Seneviratne, 2017). Here, the exceeding probability of an icing event was defined 153 as the occurrence rate (λ) where the hazard intensity s at different locations in the study area was greater than or equal 154 to the given intensity sa. That was the probability of where λ denotes the exceeding probability, sa is a given intensity of icing hazard, P(s≥sa) represents the exceeding 157 probability of the intensity that is greater than or equal to a given intensity value, and H(x, y) is the joint distribution 158 function of the indicators of inducing factors. 159 Since the hazard intensity of a disaster is determined by the different intensities of inducing factors, we can use the 160 inducing factors with different intensities to measure the hazard intensity when the disaster occurs. According to the 161 principle of factor dominance in this study, the inducing factors were solved based on the different exceeding 162 probabilities, and then the intensity values of AMP and AAT under different hazard zones were obtained. 163 Sensitivity reflects the characteristics of whether the highway is sensitive to the destructive power of a disaster, and is 176

Vulnerability of disaster-bearing bodies
proportional to the vulnerability. Given the impact of icing events on highways, this study used the grades of the 177 highway as an indicator to measure the sensitivity of the disaster-bearing body to icing events (Yang et al., 2013). In 178 terms of resistance, in areas with better economic conditions, the repair and maintenance capabilities of highways 7 the road after the disaster is higher as well. That is, the socio-economic conditions can reflect a certain extent the 181 robustness of the highways in the event of a disaster and the ability to restore to the initial state (Sugishita and 182 Asakura, 2020). Here, GDP values of different regions in the study area were used to measure the local economic 183 level. The nature of value exposure refers to the value or number of disaster-bearing bodies within the range of 184 hazards. The exposure of the social value of highways is mainly reflected in the carrying capacity of the highway per 185 unit length or area when the disaster occurs. There is no uniform standard for measuring social value. In this study, 186 congestion information, in general, was employed to measure the social value of highways. 187

Vulnerability analysis
188 Based on the above three indicators, this study described the highway vulnerability as a function of damage 189 sensitivity, resistance ability, and value exposure, which were applied to carry out the quantitative analysis of the 190 highway vulnerability (Yang et al., 2013). Here, the vulnerability of the highway was measured by the vulnerability 191 index, the greater the value, the greater the vulnerability: 192 where E is the value exposure of the highway, S represents the damage sensitivity, R is the resistance ability, Vmax and 194 Vmin represent the maximum and minimum values of the vulnerability in this study, respectively, Vs is the normalized 195 vulnerability of highways. 196 highways in the study area in the future, a risk matrix was employed here, and the highway icing risk was divided 202 into five classes using the risk score, supplemented by five colors for representation (Table 2). 203 where the red represents an extremely high risk, with the score of 1-2; Orange represents a high risk, with the score of 205 3-5; Yellow represents a medium risk, with the score of 6-10; Green represents a low risk, with the score of 12-16; 206

Risk assessment and application
Blue represents an extremely low risk, with the score of 20-25. 207 Around March 5, 12, 18, and 26, 2020, four heavy snowfalls occurred in Heilongjiang Province, causing the 208 continuous closure of many highways in the study area, including the lines of Harbin-Tongjiang, Hegang-Jiamusi, 209 and Hegang-Yichun over the study area, which seriously affected the operation of the local transportation. Taking the 210 AMP and AAT in March 2020 as an example, the hazard analysis of icing disasters was carried out, and on this basis, 211 road icing risk assessment was carried out as well. In this case, the joint distribution function of inducing factors was 212 first used to calculate the exceeding probability of icing events in the study area during this period. Subsequently, in 213 contrast to the hazard of icing events, the exceeding probability was related to the corresponding hazard zone. The The possibility of different intensities of AAT in the study area is shown in Fig. 2b The Euclidean distance between each candidate Copula function and the empirical Copula function is shown in Table  234 3. The distance between Gaussian Copula and empirical Copula is 1.2434; the distance between T-Copula and 235 empirical Copula is 0.9455; the distance between Gumbel Copula and empirical Copula is 0.0521; the distance 236 between Frank Copula and empirical Copula is 0.0892, and the distance between Clayton Copula and empirical 237 Copula is 0.1207. 238

Compared with the other four types of Copula, the Euclidean distance between Gumbel Copula and the empirical 240
Copula is the smallest, which can better characterize the dependence between AMP and AAT in the study area (Fig. 3). 241 Thus, the depiction of the dependence of inducing factors can be expressed better by the Gumbel Copula. The Gumbel density function was used here to visually describe the occurrence possibility of icing-inducing factors 248 in the study area, and further demonstrate the dependency changes (Fig. 4a) where H(x, y) represents the joint distribution function of icing-inducing factors; u, v represents the edge distribution 253 function of AMP, AAT, respectively; C represents the Copula function; x, y are the mean monthly snowfall, 254 atmosphere temperature, respectively. Indeed, the spatial surface of the joint distribution function under different 255 temperatures and precipitation is shown in Fig. 4, which contains the changing trends of Fig. 2 (a)  Here, in the principle of comprehensiveness, based on the joint distribution function of inducting factors, the hazard 262 of icing hazards is divided into five zones according to the exceeding probability (Fig. 5a). Hazard zones are from 263 high to low: Zone I, Zone II, Zone III, Zone IV, and V. 264 Based on the above zoning, the occurrence possibility, hazard zone, and intensity are linked together through the 266 exceeding probability, which intuitively presents the relationship in Fig. 5a. From the perspective of the 267 corresponding relationship, the lower the exceeding probability of a disaster, the higher the hazard intensity and the 268  (Table 4). This phenomenon can be presented intuitively 274 from the contour map of the exceeding probability (Fig. 5b). The application of the principle of factor dominance 275 can achieve the rapid assessment of icing disasters so that timely material allocation and emergency treatment can 276 be carried out for post-disaster rescue. 277 In view of the impact of icing events on highways, the sensitivity distribution of disaster-bearing bodies is measured 280 by highway grades (Fig. 6a). The higher the highway grade is, the lower the sensitivity of the road itself to icing 281 disaster. Considering the availability of data, the GDP value is used to reflect the economic conditions for a region to 282 a certain extent. GDP values of different regions in the study area are used as an indicator to measure the local 283 economic level (Fig. 6b). It is inversely proportional to vulnerability. Vulnerability is proportional to the value 284 exposure of the disaster-bearing bodies. The more congested the road, the greater the value exposure of the highway 285 when it is subjected to an icing event, and the greater the ultimate vulnerability. The value exposure of highways in 286 the study area during service is shown in Fig. 6c. 287

Damage senditivity
Grade Ⅲ Grade Ⅱ Vulnerability in this study is divided into five grades according to the natural breakpoint method in the geographical 291 information system (GIS), which can minimize the internal data of each grade and maximize the differences between 292 grades. Vulnerability grades from high to low are grade Ⅰ, grade Ⅱ, grade Ⅲ, grade IV, and Ⅴ. The higher the grade, 293 the more vulnerable it is when suffering from disasters (Table 5). 294 On the basis of grade division, GIS technology was applied to construct the vulnerability zoning map of highways to 296 explore the spatial distribution of the highway vulnerability (Fig. 7). The analysis shows that the Hegang-Yichun line 297 has the highest vulnerability in the study area. For the case application in this study, the icing hazard in the study area during this period determined by the 304 exceeding probability mainly involves three zones (Table 6). As GIS has an outperformed spatial analysis ability, it is 305 used here for the icing-hazard zoning, which can improve the accuracy of hazard assessment (Toma-Danila et al., 306 2020). Moreover, under the condition that the highway vulnerability remains unchanged, the spatial overlap about 307 Vulnerability grade Grade Ⅰ Grade Ⅱ Grade Ⅲ Grade the hazard in the case with the highway vulnerability is realized (Fig. 8a). 308 Regarding the impact of icing events on highways, after constructing the hazard distribution map of icing events, the 310 concept of the risk matrix is applied to GIS, and finally output results in the form of a risk zoning map to visually 311 display highways risk during this period (Fig. 8b) weather occurs, the interdependence between AMP and AAT will gradually increase. That is, in winter, higher 323 precipitation is generally accompanied by lower temperature, and the possibility of low temperature inducing 324 precipitation is great as well. Typically, the Gumbel Copula can better capture this information and describe it (Fig.  325   4a). 326 The risk of highway icing is related to the icing-inducing factors of different regions in the study area and the 327 distribution of highway vulnerability. Case studies have shown that as the highway vulnerability is relatively stable 328

Hazard zone
Zone Ⅱ Zone Ⅲ Zone Ⅵ Vulnerability grade Grade Ⅰ Grade Ⅱ Grade Ⅲ Grade

Risk class
Class Ⅰ Class Ⅱ Class Ⅲ Class Ⅳ Class V (a) (b) 0 250 500 125 km µ over a period of time, the icing risk encountered by highways is directly related to the possibility of icing events each 329 year, which further depends on the inducing factors of icing. In other words, we can directly assess the icing risk of 330 the highway in the future based on the forecast of temperature and precipitation in the aspect of the weather forecast. 331 This can not only complete the post-disaster deicing and relief of icing disasters but also realize the rapid assessment 332 and early warning of icing risks before the disaster. 333 It should be pointed out that the occurrence of icing disasters is not only related to precipitation and temperature but 334 also related to the wind, not considered in this study, including wind speed and direction (Somayeh et al., 2015). 335 Moreover, the quantitative analysis of disaster-bearing bodies needs to be improved in the future as well. It is not only 336 the sensitivity, resistance, and value exposure, but also the topology structure of the road network. How to embed the 337 topological vulnerability of highways into the highway vulnerability needs to be explored in the future (Faturechi and  (1) Optimal Copula determined based on the Euclidean distance showed a high upper tail dependence between 343 precipitation and temperature during icing, which could be captured by Gumbel Copula. This phenomenon showed 344 that there is a high dependence on the occurrence of low temperature and strong precipitation in extreme weather. 345 Moreover, changeable dependence could be embedded in the joint distribution function of inducing factors and 346 further used in the prediction of icing events. 347 (2) From the perspective of the law of regional differentiation, this study put forward the point of view that 348 disaster hazard zoning should be combined with the differentiation of the natural environment, which provided a 349 basis for hazard zoning. Meanwhile, taking the exceeding probability as the starting point, hazard zoning of the 350 icing was conducted, and the critical values of the hazard intensity under different zones were determined, which 351 provided a reference for the rapid assessment of icing disasters. In addition, the hazard assessment of icing disasters 352 based on exceeding probability connected the occurrence probability of icing events, the hazard zone, and intensity, 353 and formed the hazard assessment theory of icing disasters. 354 (3) The highway vulnerability constructed based on damage sensitivity, resistance, and value exposure indicated 355 that the highway with the highest vulnerability grade in the study area was the Hegang-Yichun line. The case 356 application showed that the line with the highest risk of road icing in the study area was the Hegang section of the 357 Hegang-Yichun line. Meantime, the highways around Fujin, Qitaihe, and Jiamusi also had a greater impact, which is 358 in line with the road information released by the provincial traffic center indicating the applicability of the method in 359 this study. 360 It needs to be pointed out that analysis of the highway vulnerability still needs to be improved in this study. Figure 1 Geographical location of the study area and its historical sites Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.      Vulnerability distribution of highways Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.

Figure 8
Highway icing analysis, (a) Spatial overlap of hazard and vulnerability, (b) Highway risk distribution Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.