Multi-hazard risk mapping for coupling of natural and technological hazards

Abstract Technological hazards induced by natural hazards have a significant impact on human life and economy. Their power exceeds the impact of the natural hazards that caused them. However, most of the multi-hazard risk maps failed to incorporate technological hazards when mapping multi hazard risk. To address this problem, we established a hybrid multi-hazard risk assessment methodology that calculates the risk caused by natural and technological hazards in a three-step process. First, the possible natural and technological hazards in a given region are identified by hazard-forming environment analysis and hazard source investigation respectively. A scenario simulation is adopted to analyze the intensity and impact scope of these hazards. Second, the triggered relationships among natural and technological hazards are analyzed by a hazard matrix. Then, the radar graph is adopted to calculate the multi-hazard intensity under the coupling effect of all possible hazards in this region. Third, the multi-hazard risk is calculated by aggregating land use data and multi-hazard intensity together. This methodology was applied in Huairou Science City (HSC), Beijing. The result shows that the joint action of earthquake and liquid ammonia leakage is the main reason for the high risk in the central area of HSC. The methodology developed in this study can clearly show the different impacts of each hazard on a given region under the action of hazard coupling, and the calculated results can better reflect the actual disaster situation.


Introduction
Multi-hazard risk refers to the risk arising from multiple hazards (Kappes et al. 2012;Aksha et al. 2020;Gautam et al. 2021). In principle, it accounts for the characteristics of each hazardous event (e.g. probability, frequency, magnitude), and their mutual interactions and interrelations (e.g. one hazard may occur repeatedly in time; different hazards may independently occur at the same place; different hazards may occur dependently at the same place) (Kappes et al. 2012;Marzocchi et al. 2012;Gallina et al. 2016). A multi-hazard risk map is an operational risk management tool that utilizes graphic technology to represent the identified multi-hazard risk information and visualize the development trend of risks (Kappes et al. 2012;Komendantova et al. 2014;Bathrellos et al. 2017). Since it can support the decision-making process of risk managers, it is widely applied in risk management area (Buck and Summers 2020;Pourghasemi et al. 2020).
Despite the huge amount of studies in multi-hazard risk mapping (Gallina et al. 2016;Liu et al. 2016a;Zscheischler et al. 2018), most of them merely focused on the risk induced by natural hazards, and ignored technological hazards (Barrantes 2018;Dabbeek and Silva 2020;Sutanto et al. 2020). In fact, many natural hazards have the potential to trigger technological hazards, e.g. fires, explosions and leakage of toxic or radioactive materials (Girgin et al. 2019;He and Weng 2020). These technological hazards can severely affect human life and economy, and in some cases, their impact can exceed the impact of the natural hazards that caused them (Schmidt-Thom e 2006; Chen et al. 2019;Girgin et al. 2019). For example, the occurrence of 2011 Tohoku earthquake and tsunami has triggered the meltdown of three reactors in the Fukushima Daiichi Nuclear Power Plant complex. The influence scope and duration of the nuclear leakage have exceeded those of the earthquake and tsunami (Norio et al. 2011). Therefore, in order to reduce risk, it is important to incorporate technological hazards into multi-hazard risk mapping.
Hazard interactions and interrelations analysis is the key for multi-hazard risk mapping (Tilloy et al. 2019). Previous investigations on coupling effect of natural and technological hazards mainly concentrated on the domino effect, and mostly performed event tree analysis (Marzocchi et al. 2012;Gill and Malamud 2014;Eshrati et al. 2015). Such an approach generally analyses hazard interaction that begins with given information about one primary hazard, which triggers or increases the probability of others hazards (Liu et al. 2016b). However, the interaction between different hazards is complex and dynamic, thus it is difficult to cover all situations using the linear relationship of domino effect Tilloy et al. 2019). For example, natural and technological hazards may occur independently without evident common cause, but in close proximity, spatially, temporally, or both. Hence, it is difficult to map multi-hazard risk and reflect the actual disaster situation by fully relying on the domino effect analysis. Therefore, this paper aims to develop a hybrid multi-hazard risk assessment methodology to map the risk caused by natural and technological hazards, with an explicit consideration of coupling effect between all possible natural and technological hazards.

Framework
The main aim of this paper is to develop a methodology for mapping multi-hazard risk of coupling of natural and technological hazards. The basic structure of this methodology is shown in Figure 1. The first step is to identify the possible natural and technological hazards in a given area, and analyze the intensity and impact scope of these hazards using scenario simulation method. In the second step, the overlapping areas of different hazards are identified using ArcGIS. By analyzing the possible relationships among these hazards using hazard matrix, radar graph is adopted to calculate and display the multi-hazard intensity. In the third step, land use data is selected to characterize vulnerability. A Multi-hazard risk map is drawn by combining land use data and multi-hazard intensity together.

Hazard identification and analysis
In this step, natural hazards are first identified based on the hazard-forming environment. Environmental factors in the specific geophysical environment determine the preconditions for the occurrence of a specific natural hazard (Liu et al. 2016b). According to the characteristics of these environmental factors, the spatial distribution of natural hazards in a region can be deduced. Then the historical data is adopted to simulate the possible impact scope and the corresponding intensity of these hazards, e.g. Gr€ unthal et al. (2006) simulated the inundation area and flood depth in Cologne by historic flood events. Technological hazards are identified by investigating the source of hazards, e.g. inventories of flammable and explosive materials. Then, the related models are adopted to calculate the impact scope and the corresponding intensity of technological hazards based on the hazardous material reserves, e.g. TNT equivalent method for calculating the impact scope and intensity of explosion (He and Weng 2020). Finally, the intensity of every single hazard is divided into six different intensity levels from 0 to 5. The higher the number, the greater the hazard intensity. Number 0 means that the hazard has no impact on this area.

Hazard interaction analysis
The triggered relationships among different hazards are summarized by a hazard matrix as shown in Figure 2 (Gill and Malamud 2014;Liu et al. 2016b), where hazards are categorized into natural (green) and technological (red). The vertical axis of the matrix shows the primary hazards, i.e. the initial hazards that trigger other hazards. The horizontal axis shows the potential secondary hazards, i.e. the hazards induced by initial hazards. Dark blue in the cell indicates that the primary hazard has a high probability to induce the secondary hazard. Light blue indicates that the primary hazard has a low probability of inducing secondary hazard. White indicates no obvious trigger relationship between the two hazards.  Figure 2 illustrates that most of the natural hazards have a certain probability of inducing technological hazards. 1) When an earthquake is a primary hazard, it could induce floods, tsunamis, landslides, fires, explosions, leakage of toxic materials and leakage of radioactive materials. There is a certain probability that earthquakes would result in the destruction of dams, which would cause the overflow of river and flooding (Dai et al. 2005). When an earthquake occurs on the seafloor, the abrupt rise and fall of the seabed topography would bring strong disturbance to the ocean, and there is a certain probability that it would trigger a tsunami (Norio et al. 2011). Earthquakes can easily change the stability of the slope, and it is more likely to induce landslides in a hilly area (Xie et al. 2018). There is a high probability that earthquakes would damage buildings, and cause the leakage of flammable and explosive materials, toxic materials and radioactive materials (Girgin et al. 2019). Among them, flammable and explosive materials are prone to fire or even explosion when encountering open flames (He and Weng 2020). 2) When a storm or typhoon is the primary hazard, heavy rainfall could lead to severe floods and landslides. Floods, landslides, strong winds and heavy rainfall could damage buildings, which could lead to fires, explosions and leakage of toxic or radioactive materials (Zeng et al. 2021). 3) When high temperature or drought is the primary hazard, continuous high temperature and dry climate would increase the possibility of fire and explosion (Satir et al. 2016). 4) When fire or explosion is the primary hazard, there is a certain probability that the explosion would damage the building and cause the leakage of toxic or radioactive materials (He and Weng 2020). In addition, the explosion could damage the dam and result in flooding, as well as change the stability of the slope to induce landslides.
The possible impact scope for every single hazard is identified in the hazard identification and analysis section. In this step, the overlapping area of different hazards is identified by synthesizing the impact scope of every single hazard using ArcGIS. Then, a radar graph is adopted to show the intensity degree of hazards in each assessment unit. Take Figure 3(a) as an example, the radar graph is used to show the initial intensity degree of each hazard calculated in the hazard identification and analysis section. The numbers on the graph represent the intensity degrees. With the help of radar graph, all possible hazards in a given area can be expressed and calculated. According to the triggered relationship between hazards (see Figure 2), the initial intensity degree of each hazard is adjusted before the multi-hazard intensity calculation. The adjustment principle is shown in Table 1. Using Figure 3 as an example, since earthquakes, landslides and floods all have a certain probability of inducing fires, and the intensity degree of these three hazards are 5, 2 and 3, respectively. The intensity of fire is adjusted from 2 to 4.7 (1.5 Â 1.3 Â 1.2 Â 2). Correspondingly, the intensity degrees of earthquake, landslide and flood are also adjusted as 5, 4.2, 4.3 and 6.1, respectively. The adjusted intensity degree is shown by the radar graph in Figure 3(b). Compared to the initial radar graph, the adjusted radar graph highlights the danger degree of secondary disasters (such as fires), and it can better reflect the actual disaster situation. In the hazard identification and analysis section, the initial intensity degrees of all hazards are unified to the same classification standard (0-5), and the adjusted hazard intensity degrees are also on the same classification standard. Therefore, the final multi-hazard intensity can be calculated by summing up the adjusted intensity degree of every single hazard as shown in Eq. (1). Hai (1) where, H ai represents the adjusted intensity degree of hazard i, i ¼ 1,2, … n; and MH represents the multi-hazard intensity.

Multi-hazard risk mapping
In general, the risk is defined as the probability of loss caused by the interactions between the vulnerability of exposure and the hazard (Intentional Strategy for Disaster Reduction (ISDR) 2004). Hence, in this paper, multi-hazard risk is calculated as a function of multi-hazard intensity and vulnerability of exposure (Eq. (2)). Briefly, vulnerability is given a relatively broad definition as the conditions dependent on physical, social, economic and environmental factors, which determine the potential extent of damage following exposure to hazard (Pelling 2003). Since the impact scope of a technological hazard is usually on the local scale, it is difficult to apply the For hazard that needs adjustment, there is a primary hazard that has a high probability to induce its occurrence 5 Adjusted intensity degree ¼ Initial intensity degree Â 1.5 3-4 Adjusted intensity degree ¼ Initial intensity degree Â 1.4 1-2 Adjusted intensity degree ¼ Initial intensity degree Â 1.3 For hazard that needs adjustment, there is a primary hazard that has a low probability to induce its occurrence 5 Adjusted intensity degree ¼ Initial intensity degree Â 1.3 3-4 Adjusted intensity degree ¼ Initial intensity degree Â 1.2 1-2 Adjusted intensity degree ¼ Initial intensity degree Â 1.1 vulnerability index based on administrative division statistics to the risk assessment of technological hazards. Therefore, in this paper, land use data is selected to characterize vulnerability. As shown in Table 2, according to the population density and building density of different land use types, the vulnerability is divided into five levels. For example, residential areas, schools and hospitals are the most densely populated, so residential, educational and medical land are classified as highly vulnerable areas. At last, multi-hazard risk map is drawn to show the calculated risk information with ArcGIS software.

Study area and data
Huairou Science City (HSC), covering an area of 68.4 square kilometres, is located in the northeast of Beijing. As a key national science centre in Beijing, the region contains a large number of scientific research institutions. This region is affected by the Huangzhuang-Gaoliying seismic fault zone, and the river network is relatively dense (Figure 4). In addition, there are some industrial parks in this region, which contain liquid ammonia and other hazardous materials. Therefore, the region is threatened by multiple natural and technological disasters. Three types of data are needed to implement the proposed methodology: natural hazard data, technological hazard data and land use data. Natural hazard data mainly includes meteorological data from meteorological stations, hydrologic data from hydrological stations and seismic data from China Earthquake Administration. Technological hazard data, collected from the local Emergency Management Bureau, includes the types and reserves of hazardous material. Land use data was derived from Landsat TM/ETM remote sensing image and modified through field research.

Hazard identification and analysis in HSC
According to the characteristics of hazard-forming environment (see Figure 4), e.g. the distributions of fault zones, river basins and hilly areas, the major natural hazards are earthquake, flood and landslide. Among these hazards, earthquake influences the entire study area. The closer to the fault zone, the greater the earthquake damage, thus the intensity of earthquake is indicated by the distance from the fault zone (Table 3 and Figure 5(a)) (Pacheco et al. 1993). Based on the historical rainfall data  (3) Land for municipal public facilities Medium Low (2) Forestland, cultivated land Low (1) Grassland, wasteland and hydrological data, the Hydrologic Engineering Center-River Analysis System (Brunner et al. 2016;Rawat et al. 2017) is adopted to estimate the impact scope of flood (as shown in Figure 5(b)). The inundation depth of 100-year return period rainfall is used to indicate the flood intensity (Table 3). Since soil characteristics are quite homogeneous in the study area, we only use the slope for identifying the impact scope and intensity of landslide (Table 3 and Figure 5(c)) (Guzzetti et al. 1999;Mousavi et al. 2011).
In terms of technological hazards, there are three major hazard sources in the study area, i.e. two compressed natural gas supply stations (the maximum gas storage capacity of each station is 30,000 m 3 ), and a liquid ammonia tank station (the maximum storage of liquid ammonia is 14 ton). In the case of a compressed natural gas leak, the compressed natural gas would quickly evaporate to form a vapour cloud. An explosion would occur if the vapour cloud is ignited by an open flame. Here, it is  assumed that the largest reserves of natural gas would be completely released, i.e. leakage of 30,000 m 3 natural gas in each compressed natural gas supply station. A vapour cloud explosion model is applied to estimate the damaging effect and damage radius of the explosion (Eq. (3) and (4)) (Lenoir and Davenport 1993;Zhang et al. 2019), and the intensity of explosion is classified according to the damaging effect in Table 4.
where, W TNT represents the TNT equivalent of the vapour cloud (kg); A represents the efficiency factor of the vapour cloud explosion, which is 0.04; W f is the quality of combustible gas in the vapour cloud (kg); Q f is the combustion heat of the combustible gas (kJ/kg), which is 37,000 kJ/kg; and Q TNT is the explosion heat of the TNT (kJ/kg), which is 4500 kJ/kg in this study.
where, R is the damage radius of explosion(m); and DP is the overpressure (kPa). Liquid ammonia, stored in tanks at high pressure in liquid form, evaporates into ammonia instantly after leakage. When the concentration of ammonia in the air reaches 0.5%, it would be fatal for humans to inhale for 5-10 minutes. Following previous ammonia leakage accidents (Tan et al. 2017), the release rate of leakage is set to 2.42 kg/s for 10 minutes in this study. According to the regional meteorological conditions, we set wind speed set to 0.5 m/s, 1.5 m/s, 2.9 m/s (annual average wind speed), 3.0 m/s. The atmospheric stability is set to D, E and F. Then, the Gaussian diffusion model is used to estimate the concentration distribution of the emission under different atmospheric stability and wind speed at downwind distance under the accident state (Eq. (5)) (Ji et al. 2014;Tan et al. 2017). The estimated results are shown in Table 5. The maximum distance at each concentration (bold value in Table 5) is selected to classify the intensity of ammonia leakage.
where, C(x, y, z) is the concentration (kg/s); x is the downwind distance (m); y is the transverse wind direction distance (m); z is the height above the ground level (m); Q is the release rate of leakage (kg/s); u is the average wind speed (m/s); H is the effective discharge height (m); r y and r z are the horizontal and vertical dispersion parameters, respectively.

Hazard interaction analysis in HSC
Based on the impact scope of every single hazard identified in the section 3.2, the overlapping areas of different hazards impact are identified using ArcGIS. According to the matrix in Figure 2, the possible hazard interaction types in the study area are analyzed. Earthquake has a high probability of inducing landslides, explosion, and liquid ammonia leakage, a low probability of inducing flood. Flood has a low probability of inducing explosion. Since the landslide area is far from the explosion and the liquid ammonia leakage area, it is unlikely for landslide to induce explosion and liquid ammonia leakage in the study area. The hazard interaction analysis ( Figure 6) shows that there are six hazard interaction types in the whole study area: earthquake, earthquake-flood, earthquakelandslides, earthquake-explosion, earthquake-liquid ammonia leakage and earthquakeflood-explosion. The initial radar graphs of some areas are also shown in Figure 6. Following the principle in Table 1, the radar graph is adopted to calculate the adjusted hazard intensity degree in these six types. After that, the multi-hazard intensity map is produced by synthesizing the adjusted intensity. As shown in Figure 7, the high intensity areas are mainly distributed in the core influence range of liquid ammonia leakage, and the hazard interaction type is earthquake-liquid ammonia leakage.

Multi-hazard risk mapping in HSC
According to the density of population and buildings extracted from land use data, the vulnerability index is divided into five classes (Figure 8). The final risk is calculated by aggregating the multi-hazard intensity (Figure 7) and vulnerability index (Figure 8). The final risk map is shown in Figure 9. Our analysis shows that the high risk areas are mainly distributed in the central and southwest area. The joint action of earthquake and liquid ammonia leakage is the main reason for high risk in the central area. Thus, studying the multi-hazard risk that might be caused by the coupling of natural and technological hazards is essential for risk prevention. Furthermore, the risk in the southwest area is also high due to its dense population and combined effects of earthquakes and floods.

Conclusion and discussion
The new methodology developed in this study filled a key research gap in the existing multi-hazard risk mapping research. Most of the previous multi-hazard risk maps only focused on the risk induced by natural hazards and ignore technological hazards. By contrast, the new methodology calculates multi-hazard risk with an explicit consideration of coupling effect between natural and technological hazards based on hazard matrix and radar graph. In this study, the triggered relationships among natural  and technological hazards were analyzed for the first time using a hazard matrix. The matrix can be adapted and scaled for application in specific locations to analyze the influence of the interaction between natural and technological hazards. Based on the triggered relationships shown in the matrix, the radar graph was adopted to adjust the intensity of each hazard in a specific area. The adjusted hazard intensity is the intensity under the coupling effect of all hazards in this area. In this way, the adjusted radar graph can display all possible hazards in this area, and also clearly show the impact of each hazard on the area under the action of hazard coupling.
In conclusion, with the help of the hazard matrix and radar graph, the methodology developed in this study helps public planners and decision-makers to better understand the possible interactions among different natural and technological hazards. This should allow them to implement appropriate and more targeted mitigation measures. Public planners and decision-makers can also use the assessment results to help residents and other organizations to better understand the natural or technological hazards they are exposed to, thus enhancing public risk awareness and informing local risk management.
Overall, this study has constructed a useful tool for analyzing the coupling effect between natural and technological hazards. Nevertheless, we still do not fully understand the interaction mechanism of natural and technological hazards. This study did not calculate the probability of natural hazards inducing technological hazards. Probability calculation usually requires a large amount of historical disaster data in the field of risk assessment, but historical data in most areas is insufficient. It would be interesting to explore a new way of addressing this issue without historical data in the future.

Data availability statement
Hydrologic data were collected from local hydrological stations in Huairou, and seismic data were collected from China Earthquake Administration. Technological hazard data were generated at Beijing Emergency Management Bureau. Landsat TM/ETM remote sensing image were obtained from Resource and Environmental Science and Data Center (https://www.resdc.cn/ data.aspx?DATAID=184). Derived data supporting the findings of this study are available from the corresponding author [J.F] on request.