A New Inclusive Volcanic Risk Ranking, Part 1: Methodology

The ever-increasing population living near active volcanoes highlights the need for the implementation of effective risk reduction measures to save lives and reduce the impact of volcanic unrest and eruptions. To help identify volcanic systems associated with potential high risk and prioritize risk reduction strategies, we introduce a new Volcanic Risk Ranking (VRR) methodology that integrates hazard, exposure, and vulnerability as factors that increase risk, and resilience as a factor that reduces risk. Here we present a description of the methodology using Mexican volcanoes as a case study, while a regional application to Latin American volcanoes is presented in a companion paper (Guimarães et al., submitted). With respect to existing strategies, the proposed VRR methodology expands the parameters associated with hazard and exposure and includes the analysis of 4 dimensions of vulnerability (physical, systemic, social, economic) and of resilience. In particular, we propose 41 parameters to be analyzed, including 9 hazard parameters, 9 exposure parameters, 10 vulnerability parameters and 13 resilience parameters. Since the number of parameters evaluated for each risk factor is different, they are normalized to have the same weight based on dedicated sensitivity analyses. In order to best illustrate the methodology, the proposed VRR is here applied to 13 Mexican volcanoes and compared with other approaches. We found that the volcanoes associated with the highest combination of hazard, exposure and vulnerability (3-factor VRR) for this geographic area are Tacaná and El Chichón regardless of the analyzed time window of eruption occurrence (i.e., <1 and <10 ka). Nonetheless, the volcanoes with eruption <1 ka that require the most urgent actions as associated with no or few resilience measures in place are Michoacán-Guanajuato Volcanic Field and San Martín Tuxtla (4-factor VRR); the top volcanoes in the 4-factor VRR with eruption <10 ka are Michoacán-Guanajuato Volcanic Field and Las Cumbres.


INTRODUCTION
Unplanned urbanization may lead to a concentration of population, wealth and increasingly complex economic activities in areas exposed to natural hazards, including volcanic phenomena (e.g., Zenklusen, 2007). Some highly densely populated cities are located near active volcanoes (e.g. Popocatépetl and Chichinuatzin in México, Fuji in Japan, Kelut in Indonesia, Taal in the Philippines, Vesuvius, and Campi Flegrei in Italy) increasing the associated risk; nonetheless, many of the volcanic observatories around the world still lack enough scientific and economic resources to establish resilient, multiparametric monitoring networks (e.g., Sparks 2003;Pallister et al., 2019). Because resources are often limited, monitoring and other volcanic risk mitigation strategies and efforts need to be prioritized depending on the potential impacts. Risky volcanoes can be classified based on dedicated risk ranking strategies.
Global risk rankings and maps have been proposed for various natural phenomena such as earthquakes and landslides. For instance, the Global Seismic Risk Map (v2018.1) presents the geographic distribution of average annual loss (USD) normalized by the average construction costs of the respective country (USD/m 2 ) due to ground shaking in the residential, commercial, and industrial building stock, considering contents, structural and non-structural components. The normalized metric allows a direct comparison of the risk between countries with widely different construction costs (Silva et al., 2018). In the case of landslides, a global risk map by country was constructed by Yang et al. (2015) presenting the expected annual mortality caused by these phenomena. In this approach, risk represents the interaction of hazard intensity with population density and its vulnerability. An integrative analysis of risk was presented in the World Risk Index (WRI) (Welle and Birkmann, 2015). The WRI is a statistical model for assessing the global risk of disasters that arise directly from extreme natural hazards such as earthquakes, storms, floods, droughts, or sea-level rise. It is based on the notion that disaster risk is particularly high where extreme natural events hit vulnerable societies and provides a disaster risk assessment for 181 countries worldwide. This assessment defines vulnerability as composed of susceptibility to harm, lack of coping and adaptive capacities. Regarding volcanoes, local and international efforts have been dedicated to the characterization and assessment of volcanic hazard, threat, and risk during the last decades. From the oldest to the newest, at regional scale, Bailey et al. (1983) made a general assessment of the potential for future eruptions in the United States by grouping volcanoes based on the age of the most recent eruptions and eruption periodicities. Lowenstein and Talai (1984) ranked volcanoes in Papua New Guinea, including hazard parameters based on geological features, historically recorded hazardous phenomena, and present unrest. These hazard parameters were summed with population data to calculate a relative potential hazard rating. Yokoyama et al. (1984) developed the first approach to rank volcanoes globally. This scheme assessed ten hazard parameters, focused on the characteristics of volcanism and its historical occurrences, such as magma composition, occurrence of major explosive eruption, pyroclastic density currents (PDCs), tsunamis within the last hundred years and seismic activities, and seven exposure parameters based on exposed population and historical evacuations and fatalities; these two components were scored, and the results summed to identify high threat volcanoes.
The U.S. Volcano Disaster Assistance Program developed a system of relative threat ranking for non-U.S. volcanoes (Ewert et al., 1998), which was later extended to include U.S. volcanoes by Ewert et al. (2005) and Ewert (2007) and recently updated (Ewert et al., 2018). This system considers past volcanic activity as an indication of potential future impact as well as the exposure analysis associated with each volcano. Fifteen hazard parameters were included, such as volcano type, occurrence of unrest, the general frequency of past eruptions, and the tendency toward explosivity. Unrest parameters integrate seismicity within 20 km of a volcano, deformation in response to magma intrusion or gross changes to an existing hydrothermal system and degassing. The exposure parameters include population, aviation, power generation/transmission, major development, or sensitive areas, totalizing 8 parameters. This system, with some adjustments, has been widely applied to several parts of the world, including Chile (Lara et al., 2006), Central America (Palma et al., 2009), México (Espinasa-Pereña et al., 2015, Perú (Macedo et al., 2016), Argentina (Elissondo et al., 2016) and Ecuador (Santamaría and Bernard, 2018). The threat ranking for U.S. volcanoes was recently used by Peers et al. (2021) to estimate the economic effects of volcanic alerts on the housing prices and business patterns in volcanic regions with very high-threat volcanoes.
In Japan, the Japan Meteorological Agency has developed a classification system, using two different time windows, a 10-ka activity level index, and a 0.1-ka activity level index. This system classifies volcanoes only based on volcanological criteria. In this approach exposure and vulnerability factors were not included (Uhira, 2003). To rank volcanic hazards and events that may impact the Auckland Region (New Zealand), Magill and Blong (2005a,b) calculated risk by multiplying the relative probability of the specific volcanic hazard (e.g., lava flow, lahar, ash fall) by the relative importance of the impact.
A method for volcanic hazard and exposure assessment applied to 16 priority countries of the World Bank's Global Facility for Disaster Reduction and Recovery (GFDRR) has been developed by Aspinall et al. (2011). This approach uses eight indicators to assess hazard (e.g., volcano type, Holocene PDCs, lava flow, lahar, crater lake or snow cap presence, number of sub-features, maximum Volcanic Explosivity Index -VEI, eruption frequency) and the population exposure index to evaluate the risk to the population. At hazard level, the approach attributes an uncertainty score depending on the knowledge of the volcano, which is then summed with the hazard score. Auker et al. (2015) proposed an approach (the Volcanic Hazard Index; VHI) aiming at improving the volcanic hazard assessment approaches of Ewert (2007) and Aspinall et al. (2011) and assess global volcanic hazard for the next 30 year period. This approach considers the eruption frequency and extreme characteristics (maximum VEI recorded), giving more weight to recent activity patterns, as these are considered as the most likely to indicate the type of future eruptions. Moreover, the eruption frequency is used as a multiplicator instead of just a summed indicator. In parallel, to evaluate the total volcanic threat borne by countries, Brown et al. (2015) developed two measures of volcanic threat for all volcanoes listed in the Global Volcanism Program database. First, the volcanic threat to life was obtained for each volcano based on an estimate of exposed population weighted by the number of historical fatalities. Second, the importance of volcanic threat was assessed for each volcano using an indicator based on the previous indicator multiplied by the inverse of the human development index, this index being considered as a proxy for vulnerability.
Recently, the volcanoes in Italy and Canary Islands were ranked by Scandone et al. (2016) through the Volcanic Risk Coefficient (VRC). The VRC is given by the sum of the VEI of the maximum expected eruption from the volcano, the logarithm of the eruption rate, and the logarithm of the population that may be affected by the maximum area reached by PDCs, used as a proxy to the size of the VEI. This approach is useful for comparing the degree of risk associated with different volcanoes in a region, even in case of limited availability of data. However, the lack of systematized and integrated framework still prevents us from a comprehensive and holistic understanding of the volcanic risk at regional scale.
In this paper a new strategy for Volcanic Risk Ranking (VRR) is presented that is based on the analysis of various risk factors. As all previous rankings, also the new VRR does not provide a risk assessment of volcanic systems, but a relative classification of volcanic systems based on the analysis of key risk factors. Unlike previous ranking approaches applied to volcanic settings, the new VRR not only considers hazard and exposure factors, but also physical, systemic, economic, and social vulnerabilities as well as resilience, which are the most used categories in risk analysis (e.g., Hahn et al., 2003;ISDR, 2004;Villagrán De León, 2006;Beccari, 2016;Ramli et al., 2020). In particular, volcanic systems are analyzed based on two expressions, one that combines only the factor contributing the increase the risk (i.e., hazard, exposure and vulnerability; 3-factor VRR) and one that also considers the factor that reduces risk (i.e., resilience; 4-factor VRR). The results of the two expressions are compared in order to identify volcanoes with the highest hazard, exposure and vulnerability and the volcanoes where no or only few resilience measures are in place and that, therefore, require the most urgent actions to reduce risk.
Here, we present the detailed methodology using the Mexican volcanoes as a case study. The effect of time window of analyzed volcanoes on the results is also assessed (e. g. 1 ka versus 10 ka; data provided in Supplementary Table S1). In a companion paper (Guimarães et al., submitted), the same methodology has been applied at regional scale to assess the risk of volcanoes of the entire Latin America region that have been active in the last 1,000 years.

METHODOLOGY
The first definition of volcanic risk was provided by Fournier d' Albed (1979): where risk is the possibility of a loss of life, property, or productive capacity; value is the number of people, assets (e.g. land, buildings), or economic activities (e.g. factories, power plants, agricultural land) exposed (what now is more commonly expressed as exposure); vulnerability is the expected proportion of the value be lost as a result of a volcanic event; and hazard is the probability of a specific area to be affected by a volcanic event over a certain period of time. More generally, UNDRO (1991) defined risk as "the expected number of lives lost, person injured, damage to property and disruption of economic activity due to a particular natural phenomenon": During the last few decades, several efforts have been made to integrate other key aspects into the risk analysis, such as coping capacities, which refer to the "means to face adverse consequences related to a disaster", including "management of resources before, during, and after the disaster" (Villagrán De León, 2006) or lack of preparedness (e.g., lack of early warning systems, emergency plans; Villagrán De León, 2001). Some authors integrated coping capacities at the denominator of the risk equation (e.g., ISDR, 2004;White et al., 2005;Diouf and Gaye, 2015):

Risk
Hazard × Vulnerability Coping capacity Some other authors have proposed an approach to carry out a relative comparison among systems with or without existing capacities and measures (e.g., Bollin et al., 2003;Hahn et al., 2003): where "Capacities & measures" indicate coping capacities and measures of prevention, mitigation, preparation, response and rehabilitation and reconstruction. Finally, risk has also been defined as a function of hazard, vulnerability, exposure, and resilience, without identifying a mathematical relationship between factors (MRM, 2002;Thywissen, 2006); these more qualitative analyses allow to better understand the various factors of risk, but they cannot be applied to case studies or used in a risk ranking. The new VRR integrates hazard (H), exposure (E), and vulnerability (V) as factors that contribute to increase risk (R), and resilience (Res) as a factor that reduces risk. In the proposed VRR, the resilience factor integrates both coping capacities and mitigation measures in line with the most recent definition of resilience of UNDRR (2017), i.e. "the ability of a system, community or society exposed to hazards to resist, absorb, accommodate, adapt to, transform and recover from the effects of a hazard in a timely and efficient manner, including through the preservation and restoration of its essential basic structures and functions through risk management" (https://www.undrr.org/ terminology). Therefore, resilience considers aspects that not only help to prevent and face the disasters, such as capacity to react, but also measures that can reduce the impacts, such as mitigation measures.
The selection of the equation to be used in the new VRR was based on a sensitivity analysis that tested various strategies (Supplementary Table 2A). We found that the most suited equations for a relative volcanic risk ranking are: In fact, VRR(1) (3-factor VRR) is used to investigate the factors that most contribute to increase risk, while VRR(2) (4factor VRR) is used to investigate the contribution of resilience in reducing risk. VRR(1) is in line with the risk definition of UNDRO (1991) (Eq. 2), while VRR(2) is in line with the risk definition of ISDR (2004) (Eq. 3). Nonetheless, in order to treat volcanic systems with values of resilience <1, a mathematical correction has to be made by adding the value + 1 to resilience. Eq. 6 is also preferred with respect to Eq. 4 because the multiplication allows the risk to be zero in case hazard or exposure or vulnerability is zero. In fact, based on Eq. 4, in case hazard is > 0 and exposure or vulnerability is zero, the risk would be equal to the hazard, which is not in line with the concept of risk. This can be seen, as an example, with the case of volcano Barcena, for which the hazard is > 0 and the exposure is 0 (Supplementary Table 2A); in case of Eq. 4, the risk would be equal to the hazard, while in case of Eq. 6 the risk is 0.

Hazard Parameters
Hazard is defined as a potentially damaging phenomenon that can affect a specific area over a specific recurrence time. Identification of volcanic hazards that can affect an area is typically based on records of past events, as derived from geological and/or historical and/or instrumental evidence (e.g. Alberico et al., 2002;Marzocchi et al., 2006;Connor et al., 2015;Scandone et al., 2016;Marti Molist, 2017). The hazard parameters considered in our VRR have been adapted from the methodology of Ewert et al. (2005). They concern the type of volcano, the maximum intensity, the eruptive recurrence interval and the type of volcanic phenomena. We scored the parameters related to volcano type (0 in case of cinder cone, basaltic field, small shield, or fissure vents and 1 in case of stratocone, lava dome, complex volcano, maar or caldera) to the possible type of activity such as known VEI (VEI ≥7:4, VEI 5 or 6:3, VEI 3 or 4:2, VEI 2:1 and, VEI ≤ 1:0) and to the eruptive recurrence interval between confirmed eruption, as it is reported in Global Volcanism Program (GVP). In particular, given that minor eruptions are under-recorded, the eruptive recurrence interval was analyzed only for eruptions with VEI ≥3 (eruption interval < 50 years 4; eruption interval 50-100 years 3; eruption interval 100-1000 years 2; eruption interval 1000-10,000 years or no Holocene eruption but caldera forming eruption occurred >10,000 years 1; no known Holocene eruptions 0). The occurrence of the different volcanic hazards (PDCs, lava flows, lahars, tsunamis, phreatic activity) was scored 0 (non-occurrence) to 1 (occurrence in the Holocene). The hazard associated with tephra fallout is intrinsically assessed in the VEI analysis, since the VEI scale is defined based on the volume of tephra deposits (Newhall and Self, 1982). The existence of water/ice on the volcano that could potentially trigger a lahar was also considered. The current state of activity was integrated by analyzing manifestations as seismic activity within 20 km of the volcanic edifice, since this is the common distance at which the first earthquakes appear in some observed eruptions (e.g., El Chichón, Popocatépetl, Colima; Zobin, 2012), ground deformation, and active fumarole or magmatic degassing (Table 1;  Supplementary Table S1).

Exposure Parameters
We define exposure as the identification and quantitative assessment of elements at risk (e.g., population, type of assets, infrastructure) located in an area exposed to volcanic hazards. The evaluation of the exposure of each volcano was carried out modifying and adapting the methodology of Ewert et al. (2005). For this analysis, parameters such as population, housing, critical infrastructure (e.g., transportation, power, water, telecommunication), emergency facilities (civil protection installations, police stations, fire stations, hospitals, army), critical facilities (Government offices, schools, and recreation facilities), economic activities (agriculture, livestock, forestry, fishing, mining, industry, tourism) were considered. The analysis of exposed elements for each volcano was carried out considering the possible extension of the hazards in radii of 5, 10, 30 and 100 km from the main crater. These radii were defined based on three aspects adapted from Ewert et al. (2005): 1) it was observed that in Latin America there are populations living within a radius of 5 km from the volcano crater (e.g., Tacaná volcano, in México, San Salvador volcano in El Salvador and Almolonga volcano in Guatemala), 2) the 10 and 30 km distances were chosen according to analysis in Newhall and Hoblitt (2002) for VEI 4-5 eruptions (e.g., for these magnitudes of eruption a PDC can reach 10 and 30 km distance from the main crater), 3) according to Newhall and Hoblitt (2002), tephra accumulation can exceed 10 cm at 100 km downwind in case of a VEI ≥4 eruption. A few centimeters of tephra can already have an impact on transportation, electric power, and water supply, whilst 10 cm of tephra, particularly if wet, could trigger structural damage to buildings (Ewert et al., 2005). In the case of monogenetic volcanic fields (MVF), the analysis was carried out considering all the exposed elements within the MVF as well as considering buffer distance of 5, 10 and 30 km from the rim of the field, due the variation in the spatio-temporal behavior of volcanic fields (Runge et al., 2014;Runge et al., 2015).
Population was scored based on four density intervals. Density was calculated using the ratio between the number of inhabitants and the area, considering each of the radius defined. Population density >0 inhab/km 2 , between 1 and 10 inhab/km 2 , between 10 and 100 inhab/km 2 , and above 100 inhab/km 2 score, 1, 2, 3 and 4, respectively. This density-based score system was chosen considering that the population density concept 1) simplifies the communication with stakeholders, 2) provides the opportunity to identify social conditions of a population and the pressure on the territory and on resources, and 3) allows us to fit maximum values that facilitate the normalization of the final score.
The presence of any number of residential buildings, transportation, power, water, and telecommunication infrastructures, as well as critical facilities and economic Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 697451 activities (agriculture, livestock, forestry, fishing, mining, industry, tourism) was scored with values of 1, 2, 3 and 4 considering the first appearance in 100, 30, 10 and 5 km from a central volcano, respectively or in case of MVF from a buffer (as Runge et al. (2014) and Runge et al. (2015) suggest) of 30, 10, and 5 km or first appearance within field ( Table 2).

Vulnerability Parameters
Four vulnerability dimensions were considered: physical, systemic, social, and economic, based on the most common indicators used in literature ISDR, 2004;Villagrán De León, 2006;Beccari, 2016;Ramli et al., 2020; Table 3). Physical vulnerability of buildings is defined as the expected degree of loss resulting from the impact of a certain hazard (e.g., Guillard-Gonçalves and Zêzere, 2018). Its assessment requires the evaluation of various parameters such as building category (e.g., huts, cabins, indigenous homes, or improvised dwellings), construction material, resistance and implemented protective measures. In this study, due to lack of other information, physical vulnerability analysis considers the main type of construction material (wood, masonry, and reinforced concrete) as a proxy. In particular, we assigned a score 3 in the case of wood, 2 in the case of masonry and 1 in the case of reinforced concrete. If it is not possible to know the type of material, then the typology of the houses is used as a proxy for the type of material; this criterion is based on the predominant (more than 50%) typology of residential buildings within the different radii. Most fragile buildings such as huts, ranches, cabins, mobile homes, indigenous or ethnic homes, or improvised dwellings have the highest score (3), followed by progressively less fragile structures, which score 2 such as cottage and little houses, and multiple dwelling, villa, apartment, collective dwellings that score 1. These analyses depend on available data, but it is recommended to use the latest available version of local census, whenever possible. The evaluation of the systemic vulnerability relied on evaluating how a lifeline system is prone to failure not only because of physical damage to one of its components, but also as the indirect effect of some physical, functional, or organizational failure suffered by other interconnected systems as well as on the lack of possibility to be replaced either by transfer of its function or due to redundancy. Furthermore, the method considers the vulnerability of urban and regional systems expressed as the loss of accessibility to gas, water, power, and communication utilities among others (Menoni et al., 2002). The systemic vulnerability was carried out considering two main aspects: 1) the lack of redundancy (based on the number of each type of infrastructure, e.g., transportation, power, water, telecommunication present in a  Ewert et al. (2005)  specific area) and 2) the level of accessibility (based on the number and type of transportation means available to reach an infrastructure). For the lack of redundancy, if none or only one specific infrastructure (e.g. hospital) is located within 100 km radius from the main crater or within the volcanic field up to a 30 km radius from the rim, the score is 1, otherwise it is 0 ( Table 3). The same score system was also used to assess the lack of accessibility for the same radii. Finally, considering that the management of a volcanic crisis depends on the government response, and that the presence of volcanoes on geographic borders require an effective communication between the affected countries, the existence of a volcano within 10 and 100 km from a border was scored 2 and 1, respectively. Social vulnerability refers to the inability of people, organizations, and societies to withstand adverse impacts of hazards due to characteristics inherent to social interactions, institutions, and systems of cultural values (Menoni et al., 2002). It is linked to the level of wellbeing of individuals, communities, and society. It includes aspects related to levels of literacy and education, the existence of peace and security, access to basic human rights, systems of good governance, social equity, positive traditional values, customs and ideological beliefs and overall collective organizational systems (Wisner et al., 2004). Here we consider the characteristics of population based on functional and access needs (namely proportion of children <5 years, elderly > 60 years, and invalids) as well as ethnicity, unemployment, and illiteracy (Fothergill et al., 1999;Flanagan et al., 2011;Zhou et al., 2014;Beccari, 2016;Teo et al., 2019;Federici, 2020). Such indicators are counted as the arithmetic average of data (obtained from the latest version of local census surveys) for each administrative unit (province, region, or department) inside the 100 km radius; in particular, the score value is varied from 0 to 10, with the number increasing of one unit each 10% and 10 corresponding to 100% of the population in the considered condition. The greater the proportion of  Ewert et al. (2005)).

Exposure parameters Score Maximum score value
If first appearance in 100 km radius from the main crater 1 If volcanic field, telecommunication infrastructure within the field 4 or first appearance in a buffer from the field border of 5 km 3, 10 km 2, 30 km 1 Emergency facilities: Civil protection installations, police stations, fire stations, hospitals, army If first appearance in 5 km radius from the main crater 4 4 If first appearance in 10 km radius from the main crater 3 If first appearance in 30 km radius from the main crater 2 If first appearance in 100 km radius from the main crater 1 If volcanic field, emergency facilities within the field 4 or first appearance in a buffer from the field border of 5 km 3, 10 km 2, 30 km 1 Critical facilities: Government offices, schools, recreation facilities.
If first appearance in 5 km radius from the main crater 4 4 If first appearance in 10 km radius from the main crater 3 If first appearance in 30 km radius from the main crater 2 If first appearance in 100 km radius from the main crater 1 If volcanic field, critical facilities within the field 4 or first appearance in a buffer from the field border of 5 km 3, 10 km 2, 30 km 1 Presence of at least one economic activity (agriculture, livestock, forestry, fishing, mining, industry, tourism).
If first appearance in 5 km radius from the main crater 4 4 If first appearance in 10 km radius from the main crater 3 If first appearance in 30 km radius from the main crater 2 If first appearance in 100 km radius from the main crater 1 If volcanic field, presence of at least one economic activity within the field 4 or first appearance in a buffer from the field border of 5 km 3, 10 km 2, 30 km 1 Total exposure score E 48 Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 697451  population in the considered vulnerable conditions, the greater the associated vulnerability. From an economic perspective, extreme natural events can induce damage that causes a loss of resources and impacts to a territorial system over long periods, affecting its development; such an impact can be even more significant for developing countries (e.g., Pesaro, 2018). In that respect, economic vulnerability of a country/region will depend on the different values of local resources as well as on the capability of continuing to produce goods and services and of restoring the lost resources and at which speed (Pesaro, 2018). In case economic activities are too close to the source of damage and are not If no more than one type of activity in 5 km radius from the main crater or within a volcanic field 4 4 If no more than one type of activity in 10 km radius from the main crater or within a buffer of 5 km from the border of the field in case of volcanic field 3 3 If no more than one type of activity in 30 km radius from the main crater or within a buffer of 10 km from the border of the field in case of volcanic field 2 2 If no more than one type of activity in 100 km radius from the main crater or within a buffer of 30 km from the border of the field in case of volcanic field 1 1 Total vulnerability score V 95 Monitoring provides the ability to detect and track activity frequently enough in near-real time to recognize that something anomalous is occurring. 2

Minimal
Monitoring provides the ability to detect that an eruption is occurring or that gross changes are occurring/have occurred near a volcano. Data are not collected systematically or at very long intervals (e.g., years). diversified, the economic vulnerability is considered high, affecting the capability of people and the affected region to face the consequences . If the economy is not diversified, meaning only one type of economic activity (e.g. agriculture or livestock or forestry) is present within a 5 km radius or within a volcanic field, 4 points are scored, within a 10 km radius or in a buffer of 5 km in case of volcanic field, 3 points are scored, within a 30 km radius or in a buffer of 10 km if volcanic field, 2 points are scored and within a 100 km radius or in a buffer of 30 km in case of volcanic field, 1 point is scored (Table 3).

Resilience Parameters
Resilience is defined as adaptation to changes and capacity to absorb or overcome a disturbance reaching a new level of dynamic equilibrium and/or to transform impact into opportunity (UNDRO, 1991;UNISDR, 2004;UNDRR 2017: https://www.undrr.org/terminology/resilience). It is also defined as the capacity of a system to experience disturbance and still maintain its ongoing functions and controls (Gunderson and Holling, 2002). Resilience is, therefore, considered as a positive factor that reduces the risk. Consequently, we include mitigation measures, such as the existence of hazard and risk maps, hazard-based land-use planning, monitoring, early warning systems, engineering mitigation measures, insurance coverage, educational activities for the population, exercises, and simulations both for operational institutions and for the population (Zio, 2018;Aven, 2019). Besides, we also include response capacity by considering the existence of evacuation plans, availability of shelters and, if successful, evacuations that took place in past eruptions ( Table 4). Parameters such as early warning systems or drills are score 0 (absence) or 1 (existence). Regarding hazard maps, the availability of single-hazard maps scores 1, whereas the availability of multiple-hazard maps (with different hazards described either in one single map or in multiple maps) scores 2; in fact, volcanic eruptions are typically multi hazard events (e.g., Sandri et al., 2014;Deligne et al., 2017;Takarada, 2017;Zuccaro et al., 2018;Dunant et al., 2021). The existence of risk map is given more weight and is scored 2. The score values for monitoring systems are attributed following the four categories defined by Ewert et al. (2005). Therefore, scores vary from 4 in case of a well monitored volcano, to 3 if monitoring is basic but can work in real time, to 2 if there is a limited monitoring, but still providing some key information almost in real time, to 1 if minimal monitoring is implemented, to 0 if no monitoring exists, assuming forecasting methods more accurate for the best-monitored volcanoes ( Table 4).

Risk Analysis
All the risk factors evaluated are composed of several parameters that are summed before being normalized. This normalization is based on the maximum possible value for each of the evaluated factors. In the case of hazard, the maximum value represents the highest intensity of each type of hazard; in the case of the exposure, it is the largest quantity of elements prone to be affected; in the case of vulnerability, it represents the highest level of susceptibility to damage or loss. In contrast, in the case of resilience, the maximum score represents the maximum level of capacity to face or overcome a disaster, therefore reducing the level of risk or potential impact. A total of 41 parameters were analyzed: 9 hazard parameters, 9 exposure parameters, 10 vulnerability parameters and 13 resilience parameters, where the maximum theoretical value is 19, 48, 95 and 18 for hazard, exposure, vulnerability, and resilience, respectively (Tables 1-4). The score value for each parameter was chosen after performing a sensitivity analysis with different score values (Supplementary Table 2B). Since the number of parameters evaluated for each risk factor is different, each factor was equally normalized to ten, for them to have the same weight regardless of the different scores, following standard strategies (e.g., Beccari, 2016;Ramli et al., 2020). For hazard, exposure and vulnerability, the normalization was done using the Eq. 7: Normalized factor n 1 score of all factor parameters at a given volcano n 1 maximum score of all factor parameters at a given volcano × 10 (7) where n is the number of parameters for each factor. For resilience, the normalization was achieved using the Eq. 8: Normalized factor n 1 score of all factor parameters at a given volcano + 1 n 1 maximum score of all factor parameters at a given volcano × 10 (8) It is important to notice that in Eq. 8, the mathematical factor 1 is added to resilience before normalization in order to make the equation insensitive to the normalization (see Supplementary  Table 2A).

Application of VRR(1) and VRR(2) to Volcanoes of México
VRR(1) and VRR(2) were applied to Mexican volcanoes as a case study. We considered 13 volcanoes with confirmed eruptions in the last 1,000 years, since these eruptions are the best constrained in the eruptive records (Brown et al., 2014;Mead et al., 2014). Our results were compared with those of Espinasa-Pereña et al. (2015) who applied the risk-ranking approach of Ewert et al. (2005). Finally, the influence on the results of two different time windows for eruption occurrence was also investigated considering both 1,000 and 10,000 years ( Table 5). The exposure and vulnerabilities analyses for Mexican volcanoes were conducted considering the information available in the National Atlas of Risks, an online system developed by the Civil Protection authorities that allows for the consultation of data related to exposure to different hazards (http://www.atlasnacionalderiesgos.gob.mx/ archivo/visor-capas.html) and the 2015 census of the Instituto Nacional de Geografía y Estadística (INEGI). These data were complemented by downloadable shapefile information from other institutions such as Consejo Nacional de Población (CONAPO) and from the Humanitarian Data Exchange database (HDX), the United Nations Office for the Coordination of Humanitarian Affairs (OCHA), the Humanitarian Open Street Map Team (HOT) and Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 697451 Google Street View and was processed in ArcGis 10.2. Resilience data were obtained by consulting directly with the volcano observatories. All data is provided in Supplementary Table 1.

Geological Setting and Recent Volcanic Activity
Volcanism and seismic activity in México are caused by the convergence of Cocos and North American plate in the Middle America Trench situated in the Pacific Ocean which is taking place at about 8-5 cm/year (Demets et al., 1994). In western México, volcanism is caused by the convergence of the Rivera plate beneath North American Plate about 2 cm/year (Figure 1).
In the North-western part of México there are some Pleistocene-Holocene volcanoes which origin must be related to a fossil subduction zone, but there are still debates on this issue (e.g., Negrete-Aranda and Cañón-Tapia, 2008;Calmus et al., 2011) The most recent volcanic eruption in México that resulted in catastrophic losses took place at El Chichón volcano in 1982 (VEI 5, Macías et al., 1997) which destroyed eight communities, claimed 2,000 lives, caused $117 million US Dollars of economic losses, and caused the evacuation of 22,000 people (Kreimer et al., 1999). On the other hand, Colima volcano in Western México has maintained intermittent activity since the last Plinian eruption in 1913 (VEI 4, Martin Del Pozzo and Sheridan, 1995). This condition of semi-continuous activity has dissuaded the population from urbanizing the proximal areas. The eruptions of 1997 (VEI 3) and 2016 (VEI 2) produced minimal impacts due to the existence of the monitoring system installed in the late 1980s by the University of Colima, and preventive evacuations undertaken by civil protection authorities. Popocatépetl volcano's current activity has produced inevitable damage to air navigation industry. The international airport of México City has had to cancel flights and operations during some of the most important eruptive events (Rodríguez, 2004;Guffanti et al., 2009;Nieto Torres and Martin Del Pozzo, 2017). Three preventive evacuations were carried out in association with significant activity in Additional volcanic systems could severely impact México, such as calderas and monogenetic volcanic fields (MVF). For instance, the birth of a monogenetic volcano in the Chichinautzin volcanic field represents a major volcanic threat to México City, one of the most important cultural and financial centers in Latin America and could impact the economic and social activity of the whole country (Nieto-Torres and Martin Del Pozzo, 2019).

Hazard, Exposure, Vulnerability, and Resilience
When we applied the VRR(1) and VRR(2) to Mexican volcanoes having erupted <1 ka, we found that the most hazardous volcanoes in México are Popocatépetl, Tacaná and Colima, followed by Citlaltépetl, El Chichón and Ceboruco ( Figure 2A). However, the volcanoes with the highest number of exposed elements are Jocotitlán Chichinautzin MVF, Ceboruco, Tacaná and Popocatépetl; Bárcena volcano is associated with no exposed elements ( Figure 2B). The volcanoes (with eruption <1 ka) associated with the highest vulnerability score are Tacaná, El Chichón, San Martin Tuxtla, Michoacán-Guanajuato MVF, Ceboruco and Jocotitlán ( Figure 3A). Since there are no exposed elements in Bárcena, there is also no vulnerability. Finally, of the 13 Mexican volcanoes analyzed that have had activity in the last 1 ka, only 10 are associated with some resilience parameters. Popocatépetl and Colima are the volcanoes with the highest resilience score, followed by, Tacaná, El Chichón and Chichinautzin MVF ( Figure 3B).

3-Factor VRR and 4-Factor VRR
The 4 risk factors contributing to the risk analysis are not uniformly distributed amongst the Mexican volcanoes. In particular, the controlling factors for the development of strategies to increase resilience are mostly related to the occurrence of volcanic activity within the last 1,000 years and the proximity to main cities (Figure 4). The resulting volcanoes with the highest score in the 3-factor VRR (VRR(1); Eq. 5) in México (with eruption <1 ka) are Tacaná, El Chichón, Ceboruco, Jocotitlán and Popocatépetl ( Figure 5A (Figure 5). This is mostly due to the differences in the exposure analysis between these two strategies (see Exposure for more details). The importance of resilience can be noted when the 4-factor VRR is applied (VRR(2); Eq. 6), in which case Popocatépetl, and Colima move to positions 10 and 11 because of the associated high value of resilience. In contrast, the volcanoes with the highest score in the 4factor VRR are Michoacán-Guanajuato Volcanic Field, San Martín Tuxtla and Tacaná ( Figure 5A). Tacaná is, therefore, the volcano associated with high score both in the 3-factor and in the 4-factor VRR. However, the two volcanic systems that really require urgent action to reduce risk are Michoacán-Guanajuato Volcanic Field and San Martín Tuxtla where no or only few resilience measures are in place. Similar results for the 3-factor risk ranking were obtained in case all volcanoes having erupted <10 ka were considered, with the main difference that in the 3-factor VRR (VRR(1) Eq. 5) Nevado de Toluca Volcano and La Malinche (that were not ranked in case of the <1ka window) are ranked in the first 10 positions (Nevado de Toluca Volcano is 3 rd and La Malinche is 8 th ). However, the importance of resilience can be especially noted when the 4-factor VRR is applied (VRR(2); Eq. 6). In fact, in case all volcanoes having erupted <10 ka are considered, the volcano with the highest score is still Michoacán-Guanajuato Volcanic Field but followed by Las Cumbres, La Gloria, Cofre de Perote and Valle de Bravo Volcanic Field, volcanoes with no resilience measures in place (Supplementary Figures S1, S2).

Hazard
The three most hazardous volcanoes identified based on this VRR have been permanently active since pre-Columbian times, Popocatépetl, Tacaná and Colima volcanoes. Tacaná volcano had phreatic activity in 1986. In addition, Colima volcano that ranks number 3 had a Plinian eruption VEI 4 in 1913 and has remained intermittently active to date. Moreover, its geological history shows very short periods of recurrence, only 100 years for VEI 4 eruptions, since the 16th century (Luhr and Carmichael, 1990). These volcanoes present permanent seismic activity, as well as fumaroles and, in the case of Popocatépetl and Colima, recurrent Vulcanian eruptions. Therefore, we consider these

Exposure
Jocotitlán, Chichinautzin MVF, Ceboruco, Tacaná and Popocatépetl are the volcanoes that obtained the highest exposure score. This is because density of population around these volcanoes is high, as well as the density of infrastructure that has been deployed since historical times when many populations have settled in the vicinity of these volcanoes due to the fertility of volcanic soils. Exposure analysis carried out with the VRR approach provides different results with respect to those obtained by Espinasa-Pereña et al. (2015) using the Ewert et al. (2005) approach in which Chichinautzin MVF and other monogenetic volcanic fields are on the top of the analysis. This is mainly because the Ewert et al. (2005) approach considers a logarithmic scale to quantify the number of inhabitants, unlike the VRR that considers population density. In particular, in the new VRR the number of inhabitants is analyzed based on the area in which they are distributed, which allows, particularly in the case of volcanic fields that occupy large areas, for the population not to be overestimated. Furthermore, the exposure analysis in the new VRR is carried out based on 4 different radii, while the approach of Ewert et al. (2005) considers only a radius of 30 km. In addition to the 30 km radius (from the main crater), in the new VRR we also use a 100 km radius to account for the potential impact of ashfall (in accordance with Newhall and Hoblitt (2002)), and 2 smaller radii (5 and 10 km) to better discriminate the large population densities in proximal areas that characterize volcanoes in Central America and México (Small and Naumann, 2001).

Vulnerability
Vulnerabilities are generally high in México, particularly in rural areas and in the south-eastern part of the country. The vulnerability analysis shows that the Mexican volcanoes associated with the highest vulnerability scores are Tacaná and El Chichón. These two volcanoes are in the Mexican state named Chiapas, in southeastern México, close to the border with Guatemala. In fact, Tacaná is a transboundary volcano. The State of Chiapas is the poorest state in México, a region of difficult access, where the settlements have developed as rural and farming communities lacking basic public services since historical times, with a poverty rate close to 80% (CONEVAL, 2020). This situation has generated high rates of marginalization, unemployment, inequality, as well as low educational levels and low access to information and basic services, parameters that generate vulnerability and explain the reason why volcanoes located in this territory lead the ranking of vulnerability of Mexican volcanoes.

Resilience
Popocatépetl and Colima are the volcanoes with the highest score of resilience due to the presence of many of the parameters evaluated. Popocatépetl volcano has remained active for the past  Ewert et al. (2005) approach (only based on hazard and exposure). To note that the score of Bárcena is not zero in the case of the threat assessment of Espinasa-Pereña, as in their case they had considered a certain value of exposure due to a small settlement that is no more existing.
Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 697451 26 years, whilst Colima volcano has presented intermittent activity since its last Plinian eruption in 1913. Popocatépetl, due to its proximity to México City, has become the volcano with the highest level of monitoring in México and one of the best monitored in Latin América. In addition to that, there are many other diverse activities aimed at raising awareness of the population around the Popocatépetl such as educational activities on volcanic hazards, drills, early-warning systems, design of operational plans and evacuation routes, as well as well-defined temporary shelters. These actions are sufficient for the risk score of Popocatépetl to decrease slightly and, therefore, to move from position 5 to 6 in the risk ranking ( Figure 5A). In contrast, the catastrophic eruption VEI 5 of El Chichón in 1982 and the phreatic activity of Tacaná, in 1986, have not been sufficient to implement an adequate level of monitoring in these two volcanoes and the activities of dissemination of hazard information and risk reduction have been minimal, so the level of resilience of these volcanoes is relatively low.

3-Factor VRR and 4-Factor VRR
The 3-factor VRR of Mexican volcanoes shows that Tacaná and El Chichón are the volcanoes with the highest score in México regardless of the time window considered for eruption occurrence (<1 ka or <10 ka). The high level of activity (hazard) of these volcanoes, as well as the high level of marginalization, poverty, and lack of access to the services and education present in Chiapas, as well as the low capacity for monitoring, preparedness and response make these volcanoes occupy the first levels of the risk ranking. In contrast, the case of Popocatépetl is striking, which, despite being the most hazardous volcano in México, its relatively low population density in the first 10 km radius, significantly reduces its exposure score. It also has a lower vulnerability because being in the center of México, practically all the systems analyzed have high redundancy and accessibility, populations have relatively higher socio-economic status (contrasted, for instance, with Chiapas) and have various economic activities around them. Michoacán-Guanajuato Volcanic Field, San Martín Tuxtla and Tacaná have the highest score in the 4-factor VRR. In fact, Michoacán-Guanajuato Volcanic Field has no resilience measure in place, while San Martín Tuxtla is only minimally monitored by the local university (Universidad Veracruzana) and no other resilience measures are in place. In addition, being Tacaná a border volcano, resilience measures have been put in place only on the Mexican side. Furthermore, the level of vulnerability around these three volcanoes is high; in fact, they are within the top five when considering the vulnerability analysis (Supplementary Table S1). Interestingly, Popocatépetl and Colima fall at the end of the ranking for VRR(2) because of the high score of resilience.
When the results obtained with VRR(1) and VRR(2) are compared with the results for volcanic threat carried out by Espinasa-Pereña et al. (2015), it is possible to appreciate the importance of considering the density-based score system contemplating the possible extension of the hazards in radii of 5, 10, 30 and 100 km from the crater, and other factors such as vulnerability and resilience. In fact, Popocatépetl is at the top of the threat ranking of Espinasa-Pereña et al. (2015) (similarly to our hazard results; Figure 2A). However, Popocatépetl falls to the fifth and tenth position in our 3-factor VRR and 4-factor VRR, respectively. In contrast, when vulnerability and resilience are considered, Tacaná volcano, located at the 8 th place with the Ewert et al. (2005) approach, rises to the top of both VRR(1) and VRR(2) (Figures 5A,B). This is because Tacaná has the same hazard score than Popocatépetl, but vulnerability in Tacaná is higher since it is in the poorest state of México and because the resilience associated with this volcano is low. The same situation is shown for El Chichón; this volcano ranks sixth in the hazard ranking and seventh in the exposure ranking, but second in the vulnerability ranking and low resilience which makes it occupy the second position in VRR(1) and sixth position in VRR(2) ( Figure 5A). Interestingly, Everman and Bárcena volcanoes fall in the last two positions of the ranking both for VRR(1) and VRR(2) and in the ranking of Espinasa-Pereña et al. (2015); in fact, even though they are associated with no or few resilience measures, the exposure is very low or zero. Accounting for resilience in the VRR allows to identify which volcanoes in the region have already a significant level of copying capacity and/or a series of mitigation measures in place and which, instead, require special attention for the associated risk to be decreased. Such an approach might motivate local and national stakeholders (e.g., governments) to take appropriate actions.

Influence of the Time Window of Volcanic Activity
Even though the number of Holocene volcanoes (39) is 3 times higher than the number of volcanoes that erupted during the last 1000 years (13), when these two different time windows for Mexican volcanoes were analyzed with VRR(1), no major differences in the relative ranking was observed. In fact, of the 26 added volcanoes, only Nevado de Toluca Volcano and La Malinche are ranked in the first 10 positions (Supplementary Figures S1, S2). For active volcanoes, the level of risk may not be related to the time window of volcanic activity considered but to other factors, such as the types of volcanic hazards and the number of elements exposed. Vulnerability and resilience factors are the most determining ones. Nevertheless, some differences in the ranking with VRR(2) are observed, as more volcanoes with eruption <10 ka occupy the first positions (Michoacán-Guanajuato Volcanic Field, Las Cumbres, La Gloria, Cofre de Perote and Valle de Bravo Volcanic Field). This is in line with the observation made based on Figure 4 that the occurrence of eruption within the last 1,000 years represents one of the controlling factors in developing strategies to increase resilience. As a result, many volcanoes with eruption <10 ka are associated with low resilience score and high VRR(2) score.

Mitigation Measures
The potential for major loss of life and damage to property and infrastructure of high-risk volcanoes could be reduced if key mitigation measures were implemented or, in some cases, improved. In the case of the Michoacán-Guanajuato Volcanic Field, one of the first measures that could be implemented is a formal monitoring network under the responsibility of a designated institution. Associated to this measure, a spatial-temporal volcanic Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 697451 hazard assessment should be carried out to delineate those areas within the field where an eruption is most likely to occur in the future and identify probable eruptive styles as it has been done for the Chichinautzin Volcanic Field (Nieto-Torres and Martin Del Pozzo, 2019). In the case of San Martin volcano, as the monitoring in place is minimal, a first step would be to strengthen its capability, by densifying the existing seismic network and implementing additional monitoring techniques. Moreover, it would also be important to compile both hazard and risk maps. Some volcanoes such as Tacaná, and Ceboruco, for which some mitigation measures exist, are still at the top of the 4-factor VRR for volcanoes with eruption <1 ka ( Figure 5A), as mitigation measures are not as significant as in other volcanoes (e.g., Popocatépetl). For instance, currently the Tacaná volcano, being a volcano located on the border with Guatemala, has only been instrumented from the Mexican side, while there is no monitoring network in Ceboruco. El Chichón volcano, located at the sixth place of the VRR(2), has neither a real-time seismic monitoring nor an early warning system in place and, although drills have occasionally been carried out, they do not include all the stakeholders. In addition, although hazard maps exist, they have not been translated into the languages spoken by the population, mainly indigenous.

Comparison With Previous Volcanic Ranking Approaches
Unlike previous rankings that assess mainly hazard and exposure parameters (Yokohama et al., 1984;Ewert et al., 2005;Ewert, 2007;Aspinall et al., 2011;Auker et al., 2015;Brown et al., 2015;Ewert et al., 2018) the proposed strategy was carried out by integrating 41 parameters (9 for hazard, 9 for exposure, 10 for vulnerability, and 13 for resilience) evaluated for each volcano ( Table 6). Most risk analysis methodologies use around 40 parameters (Beccari, 2016). Hazard parameters include volcano type, possible type of activity, known VEI, eruptive recurrence interval, types of volcanic hazards, and current state of activity, as was previously evaluated by Ewert et al. (2005). Some parameters evaluated in the approach of Ewert et al. (2005), such as explosive activity and major explosive activity, were excluded in our analysis since we assume that these parameters are already considered in maximum known VEI, and also because in most cases volcano eruptive history is not well known. In addition, the eruption recurrence interval was only assessed for eruptions of VEI≥3 eruptions. In the new VRR, exposure analysis consists in the integration of the density of population and housing, whereas Ewert et al. (2005) do not consider density but only consider the log number of inhabitants. Here we also consider critical infrastructure (e.g., transportation, power, water, telecommunication), emergency facilities (e.g., Civil Protection, police, fire stations, hospitals, army), critical facilities (e.g., Government offices, schools, and recreation facilities) and economic activities (e.g., agriculture, livestock, forestry, fishing, mining, industry, tourism) in a radius of 5, 10, 30 and 100 km. The distance from the volcano is also different from the previous methodologies that evaluate the exposure at 30 km from the volcano (Ewert et al., 2005).
Finally, the multiple dimensions of vulnerability and resilience had never been considered in previous rankings. Even though the number of parameters considered both for vulnerability and resilience is still limited, the new proposed strategy for both VRR(1) and VRR(2) represents an important step forward towards a comprehensive characterization of volcanic risk. In particular, resilience is still widely debated in literature (e.g., Jones, 2019;Ungar et al., 2021). Further analysis that could reconcile the various existing descriptions of resilience would largely improve the presented VRR(2).

Caveats of the New Volcanic Risk Ranking
Source and availability of data considered in the VRR can introduce significant uncertainties. The knowledge of volcanic activity and volcanic stratigraphy among volcanoes is very heterogeneous. In fact, some volcanoes have been studied and/or are monitored more than others. In addition, the potential impacts of eruptions with VEI<2 are low, and they are not typically preserved in the stratigraphic record. For these reasons, VEI≤1 scores zero while higher VEI are attributed higher scores.
The exposure and vulnerability analyses must be carried out considering as a first source of information the local census of each country and other official sources of information. When this is not possible, other free geodata could be used, such as OpenStreetMap and other similar datasets. A significant challenge is associated with these data as different types of information are not always comparable due to different scale (e.g., country level, state level or municipality level) and may also cover different time period (e.g., census). Exposure analysis is limited by the accuracy and frequency of population censuses. It is also important to consider that populations that settle in the outskirts of the cities are not always officially registered and, therefore, might not be considered in censuses; as a result, the population density analysis may underestimate the real value. Vulnerability analysis is also associated with a certain degree of uncertainty. In fact, knowing the type of material of construction and/or the type of housing in such large areas is a very timeconsuming process and the datasets used are heterogeneous, which could introduce uncertainties. Nevertheless, we consider that the criteria based on the predominant type of material and typology is adequate in case of risk ranking at regional scale. In addition, many of the lifelines are operated by private companies and many others are categorized as elements of National Security; as a result, obtaining information on their facilities is very difficult, which can cause underestimations in terms of exposure and systemic vulnerability.
Regardless of the challenges of gathering data and integrating data associated with different scales, we consider that the new VRR approach proposed here which integrates hazard, exposure, vulnerability, and resilience, represents an important tool that can be used to prioritize risk reduction strategies. It is, however, important to bear in mind that hazard, exposure, and vulnerability are dynamic factors that evolve over time at different rates. Risk analysis and ranking should, therefore, be considered as a snapshot representative of a specific time and should be updated regularly, or, at least, as soon as new data are made available. Resilience analysis is also a very dynamic factor; for instance, monitoring conditions can change very fast, when the observatories add stations, or some others are lost due to volcanic activity or because they are installed in remote places where the weather conditions are extreme. In addition, different educational activities are carried out regularly, which strengthen resilience. In some cases, it is also necessary to update the hazard and risk maps as a result the existence of a hazard and risk map does not guarantee that is representative of future activity.
As new methodologies for the characterization of risk factors and their indicators become available, they can be integrated in the VRR. This is especially true for factors such as resilience for which studies are still in their infancy. Future work could also investigate the evolution of risk factors in time, not only for the hazard as we have done in this study with a sensitivity analysis of the time window, but also for exposure, vulnerability and resilience.

CONCLUSION
The methodology of the new inclusive VRR(1) and VRR(2) allows for a relative ranking of volcanoes considering not only hazard and exposure factors, as it is the case in all previous approaches, but also the physical, systemic, economic, and social dimensions of vulnerability as well as resilience. The integration of vulnerability and resilience in the VRR provides a more accurate identification of volcanoes that require efficient risk reduction strategies.
With respect to the case study analyzed (e.g., Mexican volcanoes), Tacaná and El Chichón have scored the highest values in the 3-factor VRR (VRR(1)) regardless of the time window selected for eruption occurrence (<1 and <10 ka). This is due to their high values of hazard and exposure as well as to high values of vulnerability (physical, systemic, social, and economic). Nonetheless, Michoacán-Guanajuato Volcanic Field and San Martín Tuxtla have scored the highest values in the 4factor VRR (VRR(2)) due to the absence or only few measures of resilience, respectively. Interestingly, Tacaná also scores high in the 4-factor VRR (falling in third position) due to a combination of relatively high hazard, exposure and vulnerability and low resilience; in fact, Tacaná is the volcano with the highest score of vulnerability factor that clearly impacts both the 3-factor VRR and the 4-factor VRR ranking. Everman and Bàrcena volcanoes fall in the last two positions of both VRR(1) and VRR(2) mostly due to low or absence of exposed elements. The selection of time window of eruption occurrence (i.e. <1 and <10 ka) impacts VRR(2) more than VRR(1) as it has been observed that in México resilience measures have been mostly implemented at volcanoes with eruption within the last 1000 years.
When comparing our results with those of previous rankings for the same set of volcanoes (e.g., Ewert et al., 2005), differences are evident, particularly regarding the evaluation of exposure. The use of different radii (5, 10, 30 and 100 km) and the density-based population exposure analysis are preferred to the log of the population; in fact, in the case of extensive volcanic systems such MVF, the population parameter can be overestimated if only the amount of exposed population is considered instead of the distribution. In addition, the integration of vulnerability and resilience components is particularly important because it allows us to better identify risk reduction strategies. In fact, risk reduction strategies are most effective when based either on the reduction of vulnerability or on the increase of resilience or on both.
Since the budget for volcanic monitoring and volcanic risk mitigation is limited, especially in low-income countries, the proposed VRR methodology could be used to rank volcanic risk and, thus, help prioritize strategies and efforts, in a more efficient way to reduce the impact of volcanic activity on populations and economies. It is important to stress that this is a relative ranking that helps identify volcanoes in a region that require implementation of risk reduction strategies. Future expansion of the VRR strategy to include additional vulnerability and/or resilience parameter as they become available would improve such a prioritization.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.