Quantitative multi-hazard risk assessment to buildings in the Jiuzhaigou valley, a world natural heritage site in Western China

Abstract Mountain villages in earthquake-prone areas are likely to be exposed to multiple geological hazards. In this study, a comprehensive method is proposed to calculate the risk of buildings in two residential villages threatened by rockfalls and debris flows in the Jiuzhaigou Valley. The risk is a function of the probability of the hazard occurrence, the probability of the spatial impact, vulnerability and the economic value of the buildings. The probability of occurrence was estimated via geomorphological, geotechnical and historical record analyses; Run-out maps were obtained by simulation and used to assess the probability of spatial impact, and the vulnerability of individual buildings was calculated by both hazard intensity and building resistance. Risk curves were generated for debris flows and rockfalls for risk comparison. Uncertainties were assessed by an error propagation technique and expressed as maximum and minimum risk curves of the two hazards. Single risk maps were summed up as multi-hazard risk maps. The multi-hazard risk maps indicate that the highest risk for buildings is due to debris flows, while rockfalls pose relatively low risks to buildings. The risk result can be used as a reference for administrators and decision makers. Mitigation measures should be implemented the high-risk areas.


Introduction
Settlements with dense populations in tectonically active mountains are facing the threat of various geological hazards such as debris flows, rock avalanches in the context of climate change. The quantitative risk assessment of such multiple hazards to buildings is very important for the safety of local and resettled people, and integrated risk assessment with respect to multi-geohazards is urgently needed. The leading international organizations have realized the significance of multi-hazard risk assessment and emphasized it in their public reports (UNEP 1992;UN 2002;UNISDR 2005, UNISDR 2015. However, most of the well-established methods for the risk assessment of natural hazards have been put forward only for one kind of hazard until now, and single-hazard risk assessment is not sufficient for the elements endangered by multiple hazards. Multi-hazard risk assessment remains a challenging task due to the different characteristics of each hazard, interactive effects between different hazard processes, and diverse impacts on elements exposed in the risk Kappes et al. 2012a;Johnson et al. 2016). Many experts have conducted research on multi-hazard risk assessment in recent years. Overall, three different approaches are used to calculate compound risk (Kappes et al. 2012a). The first approach is the summation of joint indexes including hazard, vulnerability and exposure contributed by all the hazards, and the three joint indexes are commonly determined by qualitative analyses or computation indices with different weights (Dilley et al. 2005;Greiving 2006). The combination of quantitative single-hazard risk analysis method is another common approach in multi-hazard risk assessment . The two approaches have merit and demerit, the hazard interactions which can result in amplified risk or differing patterns are rarely considered (Kappes et al. 2012a). The third approach is applying cutting-edge modeling approaches (such as Bayesian networks, agent-based models, system dynamic models, event and fault trees, and hybrid models) considering the interdependencies and cascading effects (Terzi et al. 2019).
Until now, there is no universal procedure for multi-hazard risk assessment, but some quantitative multi-hazard risk approaches related to potential damage or the loss of life have been developed in some regions. Ming et al. (2015) proposed a quantitative approach of hazard risk assessment based on vulnerability surfaces and the joint return periods of hazards to assess the risk of crop losses in the Yangtze River Delta region of China. Johnson et al. (2016) carried out hazard risk assessments, including heat waves, typhoons, and landslides, in Hong Kong by utilizing indicators to describe the hazards and vulnerabilities. Chen et al. (2016) carried out a quantitative multi-hazard risk assessment in the Fella River valley, prone to debris flows and floods in the northeastern Italian Alps. The risk curve for each type of hazard was generated by the economic loss against the corresponding annual probability. In the case of earthquake-affected areas, some researchers performed preliminary studies on the multi-hazard risk assessment of various exposed elements. For example, Gr€ unthal et al. (2006) elaborated continuous building damage curves due to windstorms, floods and earthquakes in the area of Cologne (Germany) and showed the damage caused by a single hazard with different return periods.
Natural World Heritage Sites (NWHSs) are recognized as the most valuable natural asset. However, 60% of the world heritage sites are potentially exposed to at least one geological hazard such as earthquakes, landslides, volcanic eruptions, and tsunamis (Pavlova et al. 2017). Consequently, the risk assessment of geological hazards to world heritage sites is of great importance to the conservation of natural heritage (Agapiou et al. 2016;Cigna et al. 2018;Valagussa et al. 2021). Many types of research have been conducted on geological risk identification and assessment in many famous world heritage sites in recent years (Tunusluoglu and Zorlu 2009;Iriarte et al. 2010;Wang et al. 2012;Cigna et al. 2018;Valagussa et al. 2021), and many guiding principles were suggested to reduce the risks on NWHS (ICCROM UNESCO and IUCN ICOMOS 2010;Canadian Conservation Institute and ICCROM 2016). Jiuzhaigou Valley, a famous Natural World Heritage Site located in Sichuan province, southwestern China is such a typical alpine region with actively tectonic movements and highly susceptible to various geohazards that pose a serious threat to its unique landscape and safety for residents and tourists. A Ms 7.0 earthquake of with a focal depth of 20 km happened in the valley on 8 August 2017 and aggravated the situation of geological hazards. Previous studies focus on the spatial distribution of geohazards triggered by the 2017 earthquake Fan et al. 2018;Tian et al. 2019;Wang et al. 2018;Yue et al. 2018;Chen et al. 2020;Chang et al. 2021) or the risk assessment of the post-quake geohazards  Li et al. 2019). However, earthquake are high frequent in tectonically active areas. For example, Jiuzhaigou Valley is in the meizoseismal area of high magnitude earthquakes before 2017 such as the Ms 7.5 Diexi earthquake in 1933, the Ms 7.2 Songpan-Pingwu double earthquakes in 1976, and the Ms 8.0 Wenchuan earthquake in 2008. Therefore, some of the present geohazards in Jiuzhaigou are the legacy of earthquakes prior to 2017 and it needs to include previous data in a longer timespan when evaluating the multi-hazard risk. In this study, a quantitative approach is developed to assess the multi-hazard risk of  1931, 1971, 1981, 1985, 1986, 1988, 2013 and 2014 buildings in two residential villages at the Jiuzhaigou Valley More historical geohazard events are taken into account in the quantitative risk assessment to reflect the influence of historical earthquakes. The two villages, which are the main tourist accommodation and service areas, are highly susceptible to debris flows and rockfalls. The multi-hazard is likely to cause huge damage to the buildings and loss of life in the locations. Debris flow and rockfall with different scenarios are considered in the risk analysis. The risk associated with each hazard is quantified based on its intensity, spatial probability, physical vulnerability and economic value of the buildings. The intensity of hazard is modeled for these scenarios using numerical models. The study is aimed at developing a comprehensive method to evaluate the multi-hazard risk and identifying the buildings with high risk in the two villages.

Study area
Jiuzhaigou Valley is approximately 258 km south of Chengdu, the capital city of Sichuan Province and was listed in 1992 as a UNESCO World Heritage Site for its unique water landscape. It has a 43.88 km long mainstream with an average gradient of 38.7& and an elevation ranging from 1996 m above sea level (a.s.l.) in the north to 4764 m a.s.l. in the south (Figure 1). Seven deeply incised branches are distributed in the valley, with the average gradient ranging from 22% to 61.2%. The hillslope is very steep with a slope angle up to 50 in the high elevation area. The valley has a subtropical to temperate monsoon climate with a distinctive vertical variation. The annual rainfall is approximately 706 mm at Zharu village (2026 m a.s.l.), 760 mm at the Nuorilang Protection Station (2400 m a.s.l.), and approximately 1000 mm at an elevation of >3000 m. The rainfall is concentrated from May to September, mostly in the forms of heavy rains and storms. The study area is located in the transition zone from the western margin of the Sichuan Basin to the Tibetan Plateau. The exposed strata in the southern Valley is 4824.6 m-thick Upper Devonian, Middle Triassic series and in the north is mainly 2111.3 m-thick Devonian, Trias series ( Figure 2). The bedrock includes sand-slate, limestone, argillaceous limestone, and dolomite that are heavily fractured by dense faults and joints. Composite fold groups and compressed-shearing faults have been formed by intense folding and compression by many episodes of recent diastrophism (Zhang and Xu 1993). The area is located on the southern margin of the east-west of Qinling tectonic belt, east of the Songpan-Garze fold system, and adjacent to adjacent to the northeastern Longmenshan tectonic belt. The structure in the area is a complex synclinal structure with a northward tendency (Zhao et al. 2018). The Rize anticline, Jiuzhaigou fault, Zechawa fault, and Zharu fault are typical representatives of geological structures in an arcuate tectonic belt. The north-west thrust faults such as the Changhai fault, Heye fault, and Zharu fault, the north-south thrust fault such as Jiuzhaigou fault control the landform and water system of the Jiuzhaigou Valley ( Figure 2).
The Jiuzhaigou Valley lies in a seismically active region with a complex tectonic environment and scattered high slip-rate faults. According to the ground motion parameter zoning map of China, the peak ground acceleration and peak period of the seismic response spectrum in the study area are 0.20 g and 0.45 s, respectively, Since 1876, 15 earthquakes of Ms ! 6.0 have been recorded near the study area, and the magnitudes of five among them are larger than 7.0 , such as the 1933 Diexi earthquake (Ms 7.5), the 1975 Songpan triplet earthquake (Ms 7.2 and 7.2), and the 2017 Jiuzhaigou earthquake (Ms 7.0). The 2017 earthquake was caused by a concealed fault (Figure 2), triggering approximately 106 geological hazards, such as shallow landslides, rockfalls and debris flows, that pose a major threat to 17 townships and 120 villages Han et al. 2018;Lei et al. 2018;Zhao et al. 2018).
The settlement in the Jiuzhaigou Valley has a long history, and the local people are mainly Tibetan. Heye, Shuzheng Villages are two main population centers along the main tourist routes (Figures 1 and 2). The two villages are located in old debris flow deposition fans and are in serious danger of different geological hazards, such as debris flows and rockfalls, especially after the 2017 earthquake. Most of the houses in the two villages are low-rise buildings made from a variety of materials. Some buildings in the villages were partially destroyed by rockfalls triggered by the earthquake ). The multi-hazard analysis and risk assessment are fundamental for risk mitigation, especially for preventive planning, and emergency preparedness. There is an urgent need to assess the risk of buildings under multiple hazards in the villages.

Methodology
The methodology used in the study follows the main three steps: multi-hazard analysis, vulnerability assessment, and multi-risk estimation. Different types of data were collected, including digital elevation data with 5 m spatial resolution, rainfall, historical hazard events, and characteristics of the buildings in the villages. The processes of hazards with different scenarios were modeled for hazard assessment.

Hazard analysis
3.1.1. Historical hazard event Some research works have been carried out on geological hazards after the Jiuzhaigou earthquake, including field surveys, inventory mapping, and single-hazard risk assessments Chen et al. 2018;Fan et al. 2018;Hu et al. 2018;Tian et al. 2019;Ma et al. 2019). Based on these studies, more detailed information, including the type, location, and extent of damage of geological hazards around the two villages, was obtained from drone aerial photography, field surveys, and onsite interviews conducted after the earthquake (Figure 3). The two villages are threatened mainly by debris flows and rockfalls. The characteristics of each hazard are listed in Table 1.
Debris flows occur frequently and caused huge damage to the two villages in history. The 2017 earthquake triggered a great number of landslides which provide a considerable volume of loosed sediment for the debris flow ). The frequency of debris flows is likely to increase according to the statistics of post-earthquake geohazards in many areas due to largely increased sediment volume and lower rainfall thresholds (Lin et al. 2004;Cui et al. 2011;Tang et al. 2011).
Rockfalls are another type of hazard that seriously threatens the two villages. Twenty-two rockfall events were recorded from 2009 to 2013 in the Jiuzhaigou Valleys before the earthquake. The 2017 earthquake triggered 83 rockfall events, visible cracks and rock mass deformation can be observed in many hillslopes. Generally, the volumes of the arrested blocks during the earthquake are much greater than those before the earthquake. After the earthquake, the failure of the unstable rock masses still occurs due to the rainfall infiltration and their gravity. According to the investigation, two rockfall events occurred in June 2018 and August 2018, respectively. The rockfalls occurred in Heye and Shuzheng Village with small blocks prior to the earthquake, and structural countermeasures were implemented before 2017 and partially destroyed by the earthquake. The diameter of the majority of arrested blocks is between 0.02 and 0.5 m, blocks with a maximum diameter of 1.2 m crashed into a building and caused partial damage during the earthquake. After the earthquake, no rockfall events were recorded in the two villages (Institute of Mountain Hazards and Environment, CAS 2018). Based on the historical event statistics, the earthquake is the main triggering condition of the rockfalls, and the frequency of rockfall hazards is likely to increase due to the unstable rock masses near the villages.

Debris flow hazard modeling
Three debris flow events triggered by rainfalls with three different return periods (2, 10, and 50 years) corresponding to low, medium, and high magnitude are calculated in the study. The intensities of rainfall corresponding to 2-, 10-and 50-year return periods were 32, 56, and 76 mm/day, respectively. The expected volume of three debris flow hazards was calculated by the empirical approaches according to the Specification of Geological Investigation for Debris Flow Stabilization (China Geological Disaster Prevention Engineering Association 2016). The detailed calculation method for the expected volume of debris flow can be found in Liu et al. (2014).
The run-out maps corresponding to the three magnitude scenarios were obtained by numerical simulations performed using FLO-2D software (O'BRIEN 1986). FLO-2D is a simple volume conservation model that can simulate non-Newtonian flows and has been employed successfully to simulate debris flows by many researchers.

Rockfall hazard modeling
The rockfall volume is an essential parameter for hazard analysis. The estimation of the expected rockfall volume is still a challenging task. A frequency-magnitude analysis is a common method to determine the magnitude of events. However, only some records can be found along the highways in the study area (Li et al. 2019), little information about the volumes of blocks was well documented in the two villages. So we roughly estimate the block volumes and corresponding return periods according to the limit statistics of the diameters of the rock blocks falling from the slope before and during the earthquake. Blocks with three volumes for every location are considered ( Table 2). The volume of blocks with medium magnitude is equal to the maximum volume of rock blocks accumulated in the villages during the earthquake. The three magnitudes of rockfalls represent scenarios of low, medium, and high intensity, respectively. RockFall analyst was utilized to simulate rockfall propagation. RockFall Analyst is a 3D solid motion model that is capable of effectively handling large amounts of geospatial information relative to rockfall behaviors. It includes two major parts: (1) 3D rockfall trajectory simulations; and (2) raster modeling of the spatial distribution of rockfalls (Lan et al. 2007(Lan et al. , 2010. Several seeders can be uniformly sampled along a polyline according to a predefined sampling distance. DEMs, block seeders (sources), and ground surface properties are the primary inputs for calculating 3D rockfall trajectories.

Vulnerability analysis
A large number of indicators, such as building structure, number of floors, maintenance state, and surrounding condition, play a role in the potential damage of the building impacted by a hazard phenomenon. However, the importance hierarchy of the indicators varies significantly for different hazard types (Kappes et al. 2012b). Methods using indicators are still in their infancy (Papathoma-K€ ohle et al. 2017) due to hazard independence. Approaches that combine element indicators and hazard intensities have been presented by many researchers. Vulnerability models as a function of intensity (I) and resistance (R) have been proposed (Kaynia et al. 2008;Li et al. 2010;Silva and Pereira 2014). The calculation of the resistance (R) of individual buildings was based on the number of indicators mentioned above, while the intensity is related to the dynamic parameters of the hazard. In this study, the vulnerability was analyzed using the model proposed by Li et al. (2010).
where I is the hazard and R is the building resistance at risk. Both intensity and resistance range from 0 to 1. Hazard intensity is usually defined according to some parameters such as volume, velocity, and depth of the hazard. The intensity describes the damage capability to the exposed elements. For debris flow, the value of intensity was assigned based on the depth and velocity of the debris flow at the location of the specific building derived from the result of a hazard simulation (Petrascheck and Kienholz 2003;Rickenmann 2005). The classification of the intensity and proposed values are presented in Table 3. For rockfalls, the intensity was defined based on the volume of the blocks (Table 4).
Building resistance is the inherent specific characteristics of a building that can withstand a certain intensity (Li et al. 2010). All the relevant indicators should be considered (Papathoma-K€ ohle et al. 2011), but it is difficult to collect all the relevant information in practice. Related information for each building was gathered in the two villages. Four characteristics of buildings are important, including construction structure, the number of floors, building row toward the specific torrent (for debris flow) or specific slope (for rockfall), and the bounding wall range. The index importance is different and the resistance of individual buildings was calculated by using the Equation (2): where CS is the resistance scores for construction structure; NF is the resistance scores for the number of floors; BR is the resistance scores for building row toward the torrent or slope; BW is the resistance scores for the bounding wall range.     The resistance scores of the four factors of CS, NF, BR, and BW are presented in Table 5. All the values were determined based on expert knowledge by considering the historical hazards that occurred in recent years in the study area.

Economic loss estimation
Risk assessments for debris flows, and rockfall hazards have been carried out by several of researchers at different scales (Calvo and Savi 2009;Jaiswal et al. 2011;Ferlisi et al. 2012). Overall, a similar framework for risk assessment was adopted by the researchers despite the variety of hazard mechanisms and processes. In the study, the economic risk caused by a single hazard with different magnitude was summarized to obtain the total risk. Then, the economic risk caused by multi-hazard was obtained by sum the single-hazard risk. The risk assessment approach proposed by Fell et al. (2005) is used: where R(BD) is the annual damage risk of a building, P(H) is the annual probability of hazard occurrence, P(T:H) is the probability of the hazard reaching the element at risk. The P(T:H) for debris flow is determined by the run-out map and the P(T:H) for rockfalls is determined by the frequency distribution of blocks derived from the trajectory simulation. P(S:T) is the temporal spatial probability of the element at risk, V(B:S) is the vulnerability of the building at risk and E is the value of the building. When hazard events occur, the building will be definitely impacted by the hazard at the impact location because the building cannot move like people, so the P(S:T) is 1.0 in our analysis. The product of PðT : HÞ Â PðS : TÞ Â VðB : SÞ Â E is often cited as 'consequence' (Wong et al. 1997), or the product of PðHÞ Â P T : H ð Þ is often referred to as 'hazard' (Leroueil and Locat 1998). E is the economic values of the elements.
The economic value of each building was calculated using the following equation: where TA is the total area of the building (m 2 ) and MP is the reconstruction costs per square meter ($/m 2 ). The accurate MP for each building in the villages is difficult to obtain due to the various construction details and structural conditions. We estimated MP based on the field investigation of reconstruction costs in Table 6.

Debris flow modeling
Three debris flow events triggered by 2-,10-, and 50-year return period were simulated using the FLO-2D. Model parameters used in the simulation were presented in Table 7. The maximum debris flow depth and velocity during the whole process were 2600 Seeder sampling distance along the polyline 10 Seeder number at one location 10 Coefficient of normal restitution Table 9 Coefficient of tangential restitution Table 9  Friction angle  Table 9 Initial horizontal velocity (m/s) 5 Initial vertical velocity (m/s) 0 Sampling distance for seeder along polyline (m) 1 Seeder number at one location 3 derived. Figures 4 and 5 show the modeling result of the debris flow with three scenarios. The debris flow of three scenarios can inundate the majority of the buildings in Heye Village. The maximum depth in the deposition fan increases from 2.67 to 7.59 m, the maximum velocity increases from 6.35 to 13.63 m/s for the three scenarios in Heye Village. In contrast, the inundated areas of three debris flows in Shuzheng Village are quite different. The low-magnitude debris flow can only inundate a small part of the village while the debris flows with medium and high magnitude can inundate nearly the whole village. The maximum debris flow depth sharply increases from 0.33 to 4.47 m, and the maximum velocity increases from 5.59 to 10.01 m/s for the three scenarios.    The maximum debris flow depth and velocity in the location of buildings were derived to define the debris flow intensity for vulnerability calculation. The vulnerability of buildings that cannot be inundated by debris flow equals zero since debris flow poses no threat to these buildings.

Rockfall modeling
Potential release areas were based on the presence of partially detached rock blocks on the slope. Block trajectories were simulated by RockFall Analyst in the two villages. The boundary of the unstable rocks is connected as the release line for the     simulation. The simulation of rockfall propagation was performed for the unstable block along the borderline of the release area using the model parameters presented in Tables 8 and 9.
The main selected output data of rockfall consist of raster maps: trajectory with the number of blocks passing through each cell, maximum block velocity, maximum block flying height, maximum kinetic energy. Figures 6 and 7 show the results from the RockFall Analyst simulation in Heye and Shuzheng Villages. The trajectories of  large blocks with volumes of 1 m 3 on site Rockfall (1) in Heye Village are roughly consistent with the stop location of blocks triggered by the 2017 earthquake. The numbers of blocks passing through each cell reflect the probability of being struck by a falling or rolling block for buildings. The P(T:H) of buildings was calculated as the ratio of block number in the building location and the total number of blocks (Agliardi et al. 2009). The unstable rocks on site Rockfall (1) in Heye Village pose serious threat to most of the buildings with higher kinetic energy than that on site Rockfall (2).
The maximum energy for block volumes of 5 m 3 is recorded in the middle part of the slope in Rockfall (1) slope. In contrast, only a small part of buildings near the slope in the Shuzheng Villages could be struck by the blocks with high kinetic energy.

Vulnerability assessment
The characteristics of buildings were acquired from a field investigation. Buildings with resident functions are dominant in the study area. Buildings with services and commercial functions represent a large percentage and are mainly located near the main road of the area. The construction structure in the study area has distinct national characteristics in China. Figure 8 shows the statistical information of the buildings in the two villages. Figures 9 and 10 show the spatial distribution of construction structures and the number of floors of buildings in the two villages. Most buildings in Heye Village are traditional brick structures with one floor and have no bounding wall. In contrast, there are five types of construction structure in Shuzheng Village, buildings constructed with brick walls with concrete, brick, mixed brick and timber account for the majority of the total. Most buildings in Heye Village have two floors and no bounding wall. The buildings were assigned to categories and values of their resistance ability factor in Table 4. Figures 11 and 12 summarized the proportion of buildings within each resistance class. The value of single building resistance ranges from 0.2 to 1.0. Buildings with resistance to two hazards between 0.61 and 0.7 take the account greater than 50% in Heye Village. Resistance values for majority of buildings to two hazards in Shuzheng Village belong to ranges of 0.51-0.6, 0.61-0.7, and 0.81-0.9. This difference may be caused by the variety of building structures.
Building vulnerability was calculated according to Equation (1). Figures 13-16 show the building vulnerability distribution to debris flows and rockfalls with different magnitude in the two villages. Tables 10 and 11 summarized some statistics about vulnerability. The building vulnerability to debris flows in the two villages is relatively higher than that to rockfalls. For rockfall, the vulnerability of most buildings is less than 0.6. The buildings with minimum vulnerability to rockfall are mainly constructed of brick walls with reinforced concrete and consist of more than one floor with abounding wall around the whole building; the buildings are usually located relatively far away from the dangerous slope. The vulnerability to high-, medium-magnitude debris flows is greater than 0.6 for the majority of buildings in the two villages. Especially, buildings with the vulnerability of 1 to high-, medium-magnitude debris flow account for 76.80%, 60.22% of the total, respectively. These high values of vulnerability to debris flow indicate great potential damage for the buildings.

Multi-hazard risk assessment
The risks to buildings posed by a single hazard with different return periods were calculated separately, and the total risk caused by debris flow, and rockfall were calculated. The risk curves by plotting economic loss against annual probability were generated (Van Westen et al. 2014). Figure 17 displays the risk curves for economic losses, showing a large variation in risk for debris flow as compare to rockfall. The multi-hazard risk was obtained by summing up the single-hazard risks obtained by Equation 3. Figure 18 presents the economic risk maps of the two villages. The economic risk of individual buildings to the multi-hazards is summarized in Table 12. The disposable personal income (DPI) of local farmers in 2020 served as the threshold of risk classes. The average DPI of a family per year is about $10,000 in the area in 2020. When the economic damage of buildings is larger than 1000, it is a heavy loss to the family, and hence, the value was selected as the threshold of veryhigh risk. Half, 1.5 times, 2.5 times, and 3.5 times of the average DPI ($5000, $15,000, $25,000, and $35,000) is set to be the thresholds of very low, low, medium, high, and very high risks, respectively. The highest risk to a single building in the villages results from debris flow, while the risk caused by rockfall is relatively low. The maximum economic risk of a building is $55175 per year in Heye Village. A total of 18.33% of the buildings in Heye Village and 4.51% of the buildings in Shuzheng Village exhibit high or very high risk. High and very high risk-buildings in Heye Village mainly distribute in the middle and west side of the village, and about 160 residents live in these buildings. High-risk buildings in Shuzheng Village mainly distribute in the center and south area close to the main road, these buildings are mainly restaurants, shops, and hotels. Except for the residents, tourists who temporarily live in the hotels and occupy or stay in the shops are also threatened by the hazards. However, there is still a lack of distribution of tourists in the villages and the validation of the economic loss results in the two villages is difficult. Although the occurrence years of the historical events are recorded, the detailed statistics of buildings losses is still lack. In addition, most buildings in the villages were newly built in recent years and the validation of economic loss of these buildings is difficult.

Discussion
Notably, not all the buildings within the overlapping area of the multi-hazard zone (such as debris flow and rockfall) are at high risk, whereas these buildings are classified as high risk in qualitative single-hazard risk assessment . Taking into account the steps of the risk calculation, some limitations and uncertainty are involved.
The hazard frequency is determined subjectively due to the lack of long-term hazard records. The return period of debris flows is given with three scenarios in consideration of the real occurrence frequency of debris flows in history. The accurate recurrence interval of debris flows can be carried out by combined analysis of the sediment supply and rainfall threshold. The frequency distribution of rockfall volumes in the study is subjectively estimated. The lack of statistically historical catalogs of rockfall events also restricts the accuracy of the risk result. The return period estimations may be improved by dating the falling boulders, which may provide information on how many boulders of a given size have fallen in a specific period .
Concerning hazard simulation, the release areas of rockfalls are not always easy to identify and map precisely. We roughly estimated the areas on the unstable slope according to the UVA images. The presence of partially detached rock blocks and the presence of rock masses resting on unfavorably oriented joints should be used as indicators of potential release area identification. The boundary contains all the unstable blocks on the slope was used as the release lie in the simulation, which may overestimate the rockfall risk in the village. The accuracy of rockfall simulation greatly depends on the DEM cell resolution. In the study, DEM with 5 m spatial resolution is still rough for the simulation, high resolution DEM (e.g. 1 m Â 1m) may improve the quality of the simulated results.
During the vulnerability calculation process, the intensity of rockfall was assumed based on the volumes of blocks, dynamical parameters such as velocity, flying height and kinetic energy should be considered for the intensity. The uncertainty concerning the weights of indicators and resistance scores is the highest. Other factors, such as foundation and maintenance, are not included in the indicators because of data limitations. The weights of the indicators and resistance score are objectively determined. Detailed information on building characteristics and damage caused by hazards is required in future studies.
Moreover, the multi-hazard risk calculation did not consider the occurrence of the coupled events. Some buildings are simultaneously within the hazard zones of debris flow and rockfall. The economic loss of those buildings will be amplified drastically. Physical mechanical analysis of building instability is required to gain more precise results.
In order to evaluate the effect of uncertainty in the final results, a first-order second-moment procedure has been selected (Baecher and Christian 2003). The firstorder approximation of uncertainty for Equation (3) leads to: where X is the coefficient of variation, calculated as the ratio between the standard deviation (r) and the expected value (l) of the variable. The variation of P(S:T) is not considered in our analysis because the value of P(S:T) is 1.0. The X values for the different elements were performed in different ways. For the annual probability of occurrence, the vulnerability of building, and value of buildings, the X values were determined according to the limited statistics data. For the reach probability of the hazard P(S:T), the assignment of X values was based on the hazard simulation parameters. According to the statistics of historical events in the two villages, the X values of annual probability of occurrence of debris flow was assigned as 0.5, which means the return periods of the debris flow with maximum, medium, and minimum intensity vary from 25 to 75 years, 5 to 15 years, 1 to 3 years, respectively. The X values of annual probability of occurrence of rockfalls were assigned as 0.2, so the return periods of the rockfall with maximum, medium, and minimum intensity vary from 80 to 120 years, 40 to 60 years, 8 to 12 years, respectively. The X values of reach probability of the hazard were determined by the simulation results of the different hazard scenarios. The large uncertainty values were assigned to P(S:T) of debris flows, corresponding to 0.3; the uncertainty values of P(S:T) of rockfalls were assigned as 0.1. The X value of vulnerability of buildings was determined subjectively according to the intensity of hazards, corresponding to 0.5 for debris flows and 0.2 for rockfalls. The X of the value of buildings was assigned as 0.07 based on the coefficient of inflation in China. All the X values of parameters were presented in Table 13. Considering the uncertainty, the economic losses caused by debris flows in the two villages could be 61% higher or lower than the mean value shown in Figure 18, and the economic losses caused by rockfalls could be 24% higher or lower than the mean  value. Moreover, the risk curves of two hazards displaying the maximum and minimum losses according to the variation of results are present in Figures 19 and 20. The debris flow losses show a large variation between minimum and maximum losses. In contrast, the variations of rockfall losses between minimum and maximum losses are relatively small, especially in Shuzheng Village. This illustrates the large uncertainty in the debris flow risk calculation.

Conclusions
Multi-hazard risk management is indispensable for communities threatened by more than one hazard. To date, multi-hazard risk assessment remains a challenging task and does not have a standardized procedure. This study proposed a quantitative approach to calculate the risk for individual buildings exposed to the danger of debris flows and rockfalls in two densely populated villages of Jiuzhaigou Valley after the 2017 earthquake. First, the occurrence probabilities of rockfall debris flows were assumed due to the statistics of historical events. The run-out zone and related dynamic parameters of hazards were simulated in detail. Subsequently, both hazard intensity and building resistance are considered in the calculation of vulnerability. Construction structure, number of floors, building row toward the torrent or slope, and bounding wall range were adopted to determine building resistance with different weights. Finally, a quantitative risk map for a single hazard was developed, and multi-hazard risk maps for the two villages were obtained by integrating the singlehazard maps.
The final risk results show that debris flow and rockfalls pose threats with different risks to buildings in the two villages. Debris flows cause relatively high and extensive risks to buildings. In contrast, rockfalls limit the local region and only impact a small group of buildings in the villages, resulting in a relatively low risk to buildings. Overall, the high-and very-high-risk buildings are mainly located in Heye Village. The risk maps indicate that high-risk buildings are mainly located within the debris flow hazard zone with high vulnerability. Not all the buildings located in the regions affected by two hazards simultaneously belong to the high-risk class.
Despite some limitations and uncertainty in the hazard risk assessment mainly including the temporal probability of the hazard scenario, hazard modeling, vulnerability assessment and reach probability of the hazard, the multi-hazard risk map can indicate the overall risk posed to the two villages and serve as a valuable reference for the implementation of hazard countermeasures. The implementation of mitigation measures, including hazard mitigation and vulnerability reduction, is needed in the two villages, especially in the high-risk areas. Furthermore, considering a large number of tourists and residents, an assessment of risk to life should be carried out in the two villages in the future. scientific and technological support project of Land and Resources Department of Sichuan Province (Grant No. KJ-2018-24), and the CAS 'Light of West China' Program. The authors thank the two anonymous reviewers for their helpful suggestions to improve the paper.

Disclosure statement
All data included in this study are available upon request by contact with the corresponding author khhu@imde.ac.cn.