the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Simulating multihazard event sets for life cycle consequence analysis
Leandro Iannacone
Kenneth Otárola
Roberto Gentile
Carmine Galasso
In the context of natural hazard risk quantification and modeling of hazard interactions, some literature separates “Level I” (or occurrence) interactions from “Level II” (or consequence) interactions. The Level I interactions occur inherently due to the nature of the hazards, independently of the presence of physical assets. In such cases, one hazard event triggers or modifies the occurrence of another (e.g., flooding due to heavy rain, liquefaction and landslides triggered by an earthquake), thus creating a dependency between the features characterizing such hazard events. They differ from Level II interactions, which instead occur through impacts/consequences on physical assets/components and systems (e.g., accumulation of physical damage or social impacts due to earthquake sequences, landslides due to the earthquakeinduced collapse of a retaining structure). Multihazard life cycle consequence (LCCon) analysis aims to quantify the consequences (e.g., repair costs, downtime, casualty rates) throughout a system’s service life and should account for both Level I and II interactions. The available literature generally considers Level I interactions – the focus of this study – mainly defining relevant taxonomies, often qualitatively, without providing a computational framework to simulate a sequence of hazard events incorporating the identified interrelations among them. This paper addresses this gap, proposing modeling approaches associated with different types of Level I interactions. It describes a simulationbased method for generating multihazard event sets (i.e., a sequence of hazard events and associated features throughout the system’s life cycle) based on the theory of competing Poisson processes. The proposed approach incorporates the different types of interactions in a sequential Monte Carlo sampling method. The method outputs multihazard event sets that can be integrated into LCCon frameworks to quantify interacting hazard consequences. An application incorporating several hazard interactions is presented to illustrate the potential of the proposed method.
 Article
(5088 KB)  Fulltext XML
 BibTeX
 EndNote
The modeling and quantification of catastrophe risk throughout a system’s service life must encompass the occurrence of multiple natural and anthropogenic hazards. In fact, the occurrence of multiple events within a short time span (whether dictated by a causality between events or by sheer coincidence) may subject the system to exacerbated economic and societal consequences (e.g., de Ruiter et al., 2020). Such consequences have been increasing over the past decades (e.g., Di Baldassarre et al., 2018) due to several factors such as climate change, urbanization, and globalization (e.g., Cutter et al., 2008; Cutter, 2018). The assessment of such consequences cannot be approached by merely overlaying methodologies for individual hazard types. In the context of this paper, a (natural) hazard type refers to a specific category of natural events or conditions that have the potential to cause harm or damage. Indeed, multihazard risk analysis must account for the interactions among various hazard types/events and their corresponding impacts/consequences. Efforts have been devoted to standardizing the nomenclature for multihazard risk analysis (e.g., Kappes et al., 2012; Gardoni and LaFave, 2016) to enhance effective communication among end users and various stakeholders. In general, the literature refers to compound hazards (and risks) whenever there is an interaction of physical phenomena/processes (i.e., natural hazard events and their consequences) across multiple spatial and temporal scales (e.g., Pescaroli and Alexander, 2018; Sadegh et al., 2018; Zscheischler et al., 2018). However, existing research also acknowledges the importance of separating the interactions among compound hazards into occurrence interactions, which do not depend on assets/components and systems (including social ones) affected by hazard events, and impact/consequence interactions that can only occur through these exposed elements. Zaghi et al. (2016) classified the former as Level I interactions and the latter as Level II interactions. Level I interactions are due to dependencies in hazard frequencies/characteristics or the triggering or intensifying effect of one hazard type upon another. Gill and Malamud (2014) extensively reviewed Level I interactions, categorizing them based on the physical correlations between their occurrences and examining several hazard types to identify those capable of triggering or amplifying others. Nevertheless, the works mentioned above primarily categorize interactions qualitatively; they lack a discussion on the computational tools required to integrate these interactions into simulationbased frameworks for risk modeling and quantification. Consequently, the challenge of simulating sequences of events that incorporate the identified interactions largely remains unexplored. Some studies have tried to address this task. Still, they either have limited scope (e.g., sitespecific and scenariobased studies like Adachi and Ellingwood, 2008, and Marzocchi et al., 2012) or treat all interaction types uniformly, irrespective of their distinct characteristics (e.g., Mignan et al., 2014). The challenges associated with obtaining realistic sequences of events have led multiple authors to select specific, representative scenarios in their multihazard assessments, disregarding the Level I interactions in favor of a detailed study of Level II interactions (e.g., Nofal et al., 2023).
In this paper, we present a simulationbased methodology to generate sequences of hazard events (in terms of their time of occurrence and features), denoted as event sets, throughout the lifespan of a system (typically spanning 50 to 100 years, depending on the system’s socioeconomic significance). Our approach incorporates the various types of Level I interactions found in the literature, each specified by appropriate modeling techniques. We distinguish between concurrent interactions (i.e., when hazards coincide in time and space) and successive interactions (i.e., when a primary hazard precedes a secondary one). Moreover, within successive interactions, we distinguish between those where a primary hazard immediately triggers secondary events (i.e., triggering interactions) and those where a primary hazard alters the occurrence rate of secondary hazard types (i.e., altering interactions). The proposed simulationbased method assumes hazard events can be modeled as competing Poisson point processes. Nonhomogeneous Poisson processes can be incorporated by transforming them into equivalent homogeneous processes (e.g., Westcott, 1977) or leveraging thinning methods for interarrival time simulation (e.g., Lewis and Shedler, 1979), as demonstrated in this paper. The outcome is a sequential Monte Carlo (MC) approach enabling efficient simulation of multihazard event sets. These event sets can then be integrated into frameworks for Level II interactions (e.g., Selva, 2013; Dehghani et al., 2021; Otárola et al., 2024, 2023), facilitating the quantification of consequences for the purposes of life cycle consequence (LCCon) analysis.
The rest of the paper is organized as follows: Sect. 2 lists the types of hazard classifications and which dimension/information should be considered when generating multihazard event sets. Section 3 presents the proposed methodology. Section 4 shows how the procedure can be applied in practice with realistic data for an idealized casestudy location. Finally, Sect. 5 summarizes the paper’s content and proposes ideas for future work on the topic.
Hazard types can be classified based on various dimensions, influencing how they are modeled in a multihazard context. This paper specifically focuses on classifications directly influencing the simulation of multihazard event sets (i.e., a sequence of hazard events and associated features throughout the system’s life cycle) at a given target location, i.e., where there is exposure in the form of “people, infrastructure, housing, production capacities, and other tangible human assets” (UNISDR, 2005). Additional classifications could be identified based on the spatial extent, spatial variability, and spatial dependence of hazard types (e.g., Gill and Malamud, 2017). However, such classifications are outside the scope of this paper. We note that spatial considerations can be integrated into the models applied within the proposed framework (e.g., when establishing the distance between the location of interest and the location of the simulated hazard occurrence). Yet, for a more comprehensive approach to multihazard event set simulations at a larger scale (e.g., at a regional scale), one could explicitly account for various types of spatial correlations (e.g., in terms of hazard characteristics or local intensities, as defined below).
The proposed methodology employs the exceedance rates associated with different severity measures of the hazard events (see Sect. 3 for a detailed definition of a severity measure) to simulate the interarrival times between events. The proposed algorithm is agnostic toward the specific physical factors that govern the numerical values of such rates (which can be obtained from physicsbased or empirical models). Thus, we do not consider classifications like those presented in Shaluf (2007), which separate natural from humanmade hazard types. We also exclude the hazard type classifications outlined in the literature review by Gill and Malamud (2014, 2017) that refer to the specific causes of the hazard events, such as hydrological, atmospheric, or geophysical factors.
As a result, the developed algorithm incorporates the following three dimensions in the simulation of a multihazard event set (Fig. 1).

Dependency. Hazard types could be classified as independent if their occurrence and severity are not affected by the concurrent or preceding occurrence and severity of other hazard events or dependent if their occurrence and severity can be attributed to the occurrence and severity of other hazard events. In the case of dependent hazard types/events, the proposed methodology accounts for the types of interactions identified in the literature, namely concurrent and successive interactions (e.g., Zaghi et al., 2016). Concurrent interactions between two or more hazard types can be identified whenever the hazard types/events tend to occur simultaneously and/or to overlap for a period of time (e.g., storm surge, waves, and strong wind that cooccur during a hurricane). In the case of successive interactions, instead, a causal relationship exists between a primary hazard type/event and one or more secondary hazard types/events. According to these causal relationships, two broad categories can be identified within successive interactions. We denote triggering the interactions where the secondary hazard type/event (or multiple secondary hazard types/events) is triggered immediately after the occurrence of the primary hazard type (e.g., liquefaction immediately following an earthquake). In contrast, altering interactions are those where the rate of occurrence of the secondary hazard type (or multiple secondary hazard types) increases (or, more generally, changes) following the occurrence of the primary hazard type (e.g., aftershocks following a mainshock). The resulting classification of interactions is a combination of the qualitative classifications proposed by Zaghi et al. (2016) (concurrent vs. successive) and Gill and Malamud (2014) (interactions where a hazard event is triggered vs. interactions where the probability of a hazard event is increased). Section 3 describes how such classifications affect the numerical modeling of multihazard event sets.

Duration. Hazard types can be grouped into two categories based on their duration. Specifically, following the classification of disasters proposed by the Sendai Framework for Disaster Risk Reduction 2015–2030 (UNISDR, 2005), we separate suddenonset hazard types from slowonset hazard types. Suddenonset hazard types are characterized by a sudden and brief occurrence that can be modeled as a single point in time (e.g., earthquakes). Slowonset hazard types, instead, have a detectable start and end point (e.g., pandemics, droughts) and occur over an extended period.

Temporal variability. Hazard types can be classified as timeindependent or timedependent based on their rate of occurrence over time. The rate of timeindependent hazard types is constant over time (as such, the occurrence of these types of hazards is typically modeled with homogeneous processes), while the rate of timedependent hazard types varies due to physical factors that affect their probability of occurrence within a given time window (as such, the occurrence of these types of hazards is typically modeled with nonhomogeneous processes). For example, it is reasonable to assume time independence for mainshocks generating from large seismogenic zones incorporating multiple faults with similar rupture rates/characteristics (e.g., Der Kiureghian and Ang, 1977). In these cases, the occurrence of mainshocks is modeled with homogeneous Poisson processes (e.g., Abrahamson and Bommer, 2005). On the other hand, aftershocks are typically modeled with nonhomogeneous processes as their rate of occurrence typically decreases as a function of the time elapsed from the occurrence of the mainshock (e.g., Utsu, 1970). Hazard types could fall into different categories based on the complexity of the considered occurrence models. Advanced seismic hazard modeling approaches, for instance, may also consider the timevarying modeling of mainshock rates (e.g., Anagnos and Kiremidjian, 1988; Iacoletti et al., 2022), while simplified models for aftershocks may use a constant rate of occurrence that produces on average the same number of events as the nonhomogeneous process in a given time window (e.g., Iervolino et al., 2014; Iervolino and Giorgio, 2022). Timedependent hazard types can be further separated into seasonal if their rate of occurrence periodically changes over time (e.g., the occurrence of heavy rains modeled as a function of the season), increasing if their rate of occurrence increases over time (e.g., heavy rains under climate change effects), and decaying if their rate of occurrence decreases over time (e.g., aftershocks).
Figure 1 summarizes the classifications relevant to the algorithm proposed in Sect. 3. These classifications are orthogonal; i.e., they categorize hazard types based on independent criteria that do not overlap with each other. Each hazard type can be assigned to three categories based on its dependency, duration, and temporal variability.
As the primary goal of the proposed methodology is to seamlessly incorporate the available classifications for hazard interactions within a mathematical framework for hazard event simulation, the next section provides the modeling implications of the highlighted classifications.
Each hazard type is associated with certain event characteristics (e.g., rupture characteristics and magnitude in the case of earthquakes, rainfall characteristics – intensity and duration – in the case of extreme rainfall events). These quantities only characterize an event and do not account for local effects at a system's or site's location (e.g., local soil properties of a specific site, distance from the earthquake source, topography of the area). We denote the curves relating the event characteristics to their corresponding exceedance rates as event curves (e.g., magnitude–frequency distributions in seismic hazard analysis or intensity–duration–frequency in extreme rainfall event analysis). Through appropriate modeling, the event characteristics are typically translated into site and/or systemspecific intensity measures (e.g., ground shaking for structures with selected vibration periods, flood depths), which are typically used as an input to obtain the corresponding physical impacts caused by the event at a given location (e.g., through fragility/vulnerability models; e.g., Gentile et al., 2022). Appropriate methods for the local intensity calculation can be found in the literature and depend on the hazard type considered. For example, ground motion models (GMMs) (e.g., Douglas and Edwards, 2016) can be used to translate earthquake characteristics into earthquakeinduced ground motion intensity measures such as peak ground parameters (i.e., peak ground acceleration, velocity, and displacement) and pseudospectral accelerations, among others. Such models could also account for the spatial and crossintensity correlation of the intensity measure (e.g., Jayaram and Baker, 2010a, b). For floods, accurate flowbased hydraulic models can be used to translate the rainfall characteristics into flood depths at different locations (e.g., Mignot and Dewals, 2022). In general, analyses at the regional scale call for maps displaying the spatial variability in the intensity measure across different locations. Finally, the methods to generate intensity measures can be integrated within endtoend probabilistic frameworks to obtain curves linking each intensity measure value to its associated exceedance rate in a given time window (i.e., 1 year in the case of annual exceedance rates). Such curves are denoted as hazard curves (e.g., the curves for the exceedance of a given wind speed and surge depth in Apivatanagul et al., 2011, or the hazard curves from the Global Earthquake Model Seismic Hazard Map; Pagani et al., 2020). Because both event and hazard curves are used interchangeably in the proposed formulation for the same purpose on a casebycase basis, we arbitrarily introduce the term rate curves to refer collectively to both cases and the term severity measure to refer to both event characteristics and intensity measures. In fact, some aspects of multihazard event sets may be governed by the event characteristics (i.e., the rate of aftershocks is governed by the magnitude of the mainshock). In contrast, others may be governed by the intensity measures (i.e., the triggering of a landslide following the occurrence of an earthquake is governed by the ground motion at the slope location).
3.1 Mathematical modeling of event sets
Let us focus first on independent, suddenonset, timeindependent hazards. The discussion will then be extended to incorporate the dependencies across events and to account for slowonset and timedependent hazards. The severity measure associated with the occurrence of the ith hazard type h_{i} is denoted as m_{i}, and the corresponding mean exceedance rate (exceedance rate for brevity hereinafter) is denoted as λ(m_{i}). It must be stressed that, as the hazard types are timeindependent, such rates do not change as a function of time. The mean occurrence rate (occurrence rate for brevity hereinafter) of hazard type h_{i} can be obtained from the rate curve as
where m_{i,min} is the minimum value of interest of the severity measure (e.g., for earthquakes, it could be the minimum magnitude of engineering interest). In other words, the occurrence rate of h_{i} is the exceedance rate of its minimum severity measure. A schematic representation of a rate curve can be seen in Fig. 2b. If a hazard type is associated with multiple severity measures, rate surfaces define their joint exceedance rate. For example, intensity–duration–frequency surfaces are a standard tool to quantify the mean return period of given rainfall heights and durations (e.g., Fadhel et al., 2017). They define the mean return period (reciprocal to the rate in the case of timeindependent hazard types) as a function of both severity measures.
Slowonset events are defined by the rate curve associated with the start of the event, ${\mathit{\lambda}}_{i}^{\mathrm{s}}\left({m}_{i}\right)$, and the rate associated with the end of the event, ${\mathit{\lambda}}_{i}^{\mathrm{e}}$. The rate curves for timedependent hazard types are instead described as a function of time.
The rates obtained from these curves/surfaces are used in event simulation, assuming that the event occurrences follow a Poisson process, either homogeneous or nonhomogeneous for timeindependent and timedependent hazard types, respectively (for example, the Bartlett–Lewis and the Neyman–Scott models for storm generation in Ritschel et al., 2017). For homogeneous Poisson processes, the interarrival times t_{h} between event occurrences of hazard type h_{i} follow an exponential distribution with parameter λ_{i}. In this case, simulating a multihazard event set consists of randomly sampling numbers from an exponential distribution. A critical assumption of homogeneous processes is that events occur independently of each other, a somewhat restrictive assumption for specific hazard types affected by seasonality and/or previous hazard occurrences (the hazard types classified as dependent and timedependent in Sect. 2). To account for such hazard types, we simulate events from nonhomogeneous Poisson processes using a procedure known in the literature as thinning (Lewis and Shedler, 1979). This procedure (described in detail in Sect. 3.3.2) also relies on the random sampling of exponentially distributed numbers and can be efficiently implemented in practice.
3.2 Required input
Figure 2a shows a portion of the interaction matrix from Zaghi et al. (2016), which includes flood (F), heavy rain (HR), earthquake (E), and landslide (L). The classification between concurrent and successive interactions is kept unaltered from the original reference. However, the successive interactions have been further separated as triggering (L → F, F → L, HR → L, E → L) and altering (E → E, L → L). The distinction between triggering and altering is needed to capture the different implications of these interactions in the modeling framework.
Each type of interaction requires different information to be modeled, provided in the list below and shown in Fig. 2. For successive interactions, we provide descriptions for the case when h_{2} is a secondary hazard type following the occurrence of a primary hazard type h_{1}. Depending on the considered interaction, the same hazard type could be classified as primary or secondary (for example, for the set of hazard types in Fig, 2, an earthquake could be a mainshock – primary – or an aftershock – secondary). Also, the discussion can be extended to the case with multiple secondary hazard types.

Concurrent interactions (Fig. 2d) are defined by the joint rate of exceedance of the severity measures of all hazard types involved (e.g., the joint exceedance of a given snow depth and a given wind speed, as in Wang and Rosowsky, 2013). This results in rate surfaces that can be interpreted analogously to the ones for single hazard types described by multiple severity measures.

Successive triggering (Fig. 2c) interactions are defined by the probability of occurrence of h_{2} conditioned on the severity of h_{1}, i.e., P(h_{2}m_{1}). In cases where the severity measure of h_{2}(m_{2}) is of interest, a conditional probability distribution of such quantity (f(m_{2}m_{1})) is also provided (conditioned on the severity measure of h_{1}). As the secondary hazard event(s) is assumed to occur immediately or shortly after the occurrence of the primary hazard event, there is no time component incorporated into the modeling of triggering interactions. An example of conditional probabilities and conditional distributions used to model triggering interactions can be found in Neri et al. (2008), where the probability of several secondary hazard events, such as floods and landslides, is conditioned on the occurrence of the volcanic eruption of Mount Vesuvius. Neri et al. (2008) also provide the variability in the severity measures associated with secondary hazard events. Another example can be found in Parker et al. (2015), where the authors quantify the probability of a landslide on a slope conditioned on the occurrence of an earthquake and its severity.

Successive altering (Fig. 2e) interactions are defined by the change in the rate curves of h_{2} following the occurrence of h_{1}. Mainshocks and aftershocks are an example of successive altering interactions. Following a mainshock, the rate of aftershocks is typically modified in terms of the characteristics of the mainshock using the modified Omori law (Utsu, 1970) or more advanced models (e.g., Iacoletti et al., 2022). Rates for this type of interaction typically decay with time as the memory of the primary hazard subsides, resulting in nonhomogenous Poisson processes. This paper proposes a thinning methodology to incorporate nonhomogeneous processes into the formulation (Sect. 3.3.2). Alternatively, the nonhomogeneous processes can be translated into an equivalent homogeneous process with a given memory mem_{i}. By taking the inverse of such memory, we can also define a “memory loss rate”, ${\mathit{\zeta}}_{i}=\mathrm{1}/{\mathrm{mem}}_{i}$. This rate determines how often the system loses its memory of h_{1} and the occurrence rates of h_{2} return to the original level (e.g., zero for aftershock occurrences). We call rate curves unaffected by altering interactions original rate curves and rate curves affected by altering interactions conditional rate curves.
3.3 Life cycle hazard event set simulation
We incorporate the above information into a sequential MC simulation approach. The procedure outputs life cycle hazard event sets, i.e., multihazard event sets (times of occurrence, hazard types, and associated severity measures), throughout the life cycle of the system of interest. Because of the simplifying assumptions, the simulation of such event sets is computationally efficient while retaining the relevant implications of the hazard interactions. It can be repeated multiple times to obtain relevant statistical quantities such as (i) the probability of having a given number of occurrences of a specific hazard type in a given time span; (ii) the probability of occurrence of a specific combination of hazard events within the life cycle; and (iii) the probability distribution of the severity measures of the hazard events and joint probability distributions of the severity measures of multiple hazard events occurring within a short time frame. These quantities, as well as the simulated event sets, can be incorporated into formulations for Level II interactions (e.g., Selva, 2013; Dehghani et al., 2021; Otárola et al., 2024, 2023) to obtain the expected consequences of the hazard events throughout the life cycle of a system. To the best of the authors' knowledge, no algorithm is currently available in the literature that accounts for the types of interactions and additional aspects (i.e., event duration and temporal variability) highlighted in this paper. An alternative sequential MC approach has been proposed by Mignan et al. (2014), which separates the simulation of primary events from the simulation of secondary events (all primary events are simulated, and then all secondary events are simulated). However, such an algorithm only considers suddenonset, timeindependent hazard events and models all interactions as successive triggering. Similarly, Selva (2013) used simplified, closedform solutions to translate rate curves into probabilities of occurrence of the hazards within a given time period. While the interactions between hazards can be included by modifying the probability of occurrence of the secondary hazard (Selva, 2013, introduces “coactive risk factors” for this purpose and provides an example with volcanic eruptions and ash fallout), such an approach is also limited to suddenonset, timeindependent events and does not capture the intricacies of hazard sequences that may include multiple successive interactions.
The following subsections describe the proposed algorithm to generate the event sets, starting with the case with timeindependent, suddenonset hazard types and then extending the discussion to timedependent and slowonset hazard types. It is worth noting that, in the presence of a single hazard (e.g., earthquakes) from different sources, the occurrence of the event from each source can be modeled as its own hazard type (e.g., h_{1} is earthquake from source 1, h_{2} is earthquake from source 2) each characterized by a given rate of occurrence (and characteristics).
3.3.1 Proposed algorithm
Figure 3 shows a visual representation of the proposed algorithm to simulate event sets for the case with timeindependent, suddenonset hazard types. The flowchart in Fig. 4 details the described sequential MC approach, with reference to each step shown in Fig. 3. Every type of interaction is included in the proposed procedure and incorporated based on its specific characteristics.
The simulation of the event set starts in a neutral state where the rates for each hazard type ${h}_{i}(i=\mathrm{1},\mathrm{\dots},N)$ are defined based on the corresponding original rate curves (Step 1 in Fig. 3). We define I as the set of indices of the hazard types that, at any given time, have affected the rate of any secondary hazard type (e.g., if the system retains memory of hazards h_{1} and h_{3}, we will have I≡{1,3}). Because the system has no memory of any previous hazard events at this point, we set I=∅. The theory of competing Poisson processes determines that the rate of occurrence of the first hazard event is equal to the sum of the rates of the individual hazard types. Consequently, the time of occurrence t of the first hazard event is sampled from an exponential distribution (${f}_{T}\left(t\right)=\mathit{\lambda}{e}^{\mathit{\lambda}t},t>\mathrm{0}$) with parameter $\mathit{\lambda}={\sum}_{n=\mathrm{1}}^{N}{\mathit{\lambda}}_{n}$ (Step 2 in Fig. 3). Once the event has been simulated, it is assigned to one of the ith hazard types h_{i} (Step 3 in Fig. 3). The probability that the hazard event belongs to the ith hazard type h_{i} is
Three phases follow the simulation of a hazard event: Phase 1 is the assessment of the hazard severity (i.e., m_{i}), Phase 2 is the simulation of hazard events caused by triggering interactions, and Phase 3 is the reassessment of the rates based on altering interactions.
In Phase 1, the severity measure of the simulated hazard event m_{i} is obtained from the rate curve of the ith hazard type λ_{i}(m_{i}) using the onetoone relationship between the cumulative distribution function (CDF) of the random variable M_{i} and its rate of exceedance, as
In the case of hazard types associated with multiple severity measures and/or concurrent hazards, all severity measures are obtained from rate surfaces rather than rate curves. A similar relationship to Eq. (3) can be obtained for the multidimensional case. Phase 1 is not shown in Fig. 3 for visual clarity.
In Phase 2, the triggering interactions are simulated (Step 4 in Fig. 3). The occurrence of each secondary hazard event is simulated by considering the conditional probability established within the triggering interaction’s definition. Subsequently, the severity measure of the secondary event(s) is sampled from the corresponding conditional probability distribution.
In Phase 3, the rates of each hazard type are reassessed to account for the altering interactions. In particular, for each altering interaction (Step 5 in Fig. 3) (i) we substitute the original rate curves for each secondary hazard type with the corresponding conditional rate curve; (ii) we introduce an additional “memory loss” Poisson event with rate ${\mathit{\zeta}}_{i}=\mathrm{1}/{\mathrm{mem}}_{i}$ (associated with the primary hazard event) to the pool of possible events; and (iii) if i∉I, we add i to the set of indices I.
We can then simulate the following event, which can be the occurrence of either a new hazard event (with rate λ_{i}) or a memory loss event (with rate ζ_{i}). The theory of competing Poisson processes determines that the rate of the next event occurrence is equal to the sum of the rates of the individual events (which now include both hazard events and memory loss events). Consequently, the time of occurrence of the next event is sampled from an exponential distribution with parameter $\mathit{\lambda}={\sum}_{n=\mathrm{1}}^{N}[{\mathit{\lambda}}_{n}+{\mathit{\zeta}}_{n}{\mathbf{1}}_{\mathit{\left\{}n\mathrm{\in}\mathit{I}\mathit{\right\}}}]$ (Step 6 in Fig. 3), where 1_{{∙}} is the indicator function (i.e., 1_{{nϵI}}=1 if n∈I and 1_{{nϵI}}=0 if n∉I). The simulated event is then assigned either to one of the hazard types ({H=h_{i}}) or to the loss of memory of one of the hazard types ($\mathit{\{}H={h}_{i}^{*}\mathit{\}},i\mathit{\u03f5}\mathit{I}$) (Step 7 in Fig. 3). The probability that the next event is the occurrence of the ith hazard type ({H=h_{i}}) is
and the probability that the next event to occur is a loss of memory of the ith hazard type ($\mathit{\{}H={h}_{i}^{*}\mathit{\}},i\mathrm{\in}\mathit{I}$) is
If the simulated event is a hazard event, Phases 1–3 are repeated. If the simulated event is the memory loss of the ith hazard type, we remove i from the set of indices I, replace the corresponding conditional rate curves with the original rate curves, and remove the Poisson event with rate ζ_{i} from the pool of possible events (Step 8 in Fig. 3).
3.3.2 Incorporating timedependent events
This section details how to modify the procedure in Sect. 3.3.1 to incorporate nonhomogeneous Poisson events used to capture the temporal variability in the hazard (timedependent events in Sect. 2). We focus on the case with multiple decaying processes associated with altering interactions (e.g., aftershocks after the occurrence of mainshocks). The procedure can be adapted to the case of seasonal processes with slight modifications. Figure 5 visualizes the described procedure. To aid in the interpretation of Fig. 5, the reader may think of hazard type 1 as mainshocks and hazard type 3 as aftershocks.
In this case, we define J as the set of indices of the (secondary) hazard types that have been affected by altering interactions. No memory loss events need to be considered in this case. Because all the considered nonhomogeneous processes follow from altering interactions and the original rate curves are associated with homogeneous processes, Steps 1–4 described in Sect. 3.3.1 are unchanged. After the occurrence of a primary hazard event associated with at least one altering interaction, we define the timevarying rate of the nonhomogeneous process associated with the ith (secondary) hazard type as λ_{i}(t−t_{0}) (for each of the secondary hazards of the altering interactions), where t_{0} is the time of occurrence of the primary hazard event that caused the altering interaction (Step 5 in Fig. 5). We also add the indices associated with the secondary hazard types to J, if they were not already included in the set (i.e., i∉J). To simulate the occurrence of the next hazard event, which results from multiple competing homogeneous and nonhomogeneous processes, we use a modified, sequential version of a procedure known in the literature as thinning, which involves using a higher, homogeneous rate for event simulation and then discarding a selection of the simulated events, i.e., classifying them as null events (Lewis and Shedler, 1979). In the proposed sequential approach, for each i in J, we fix the rate of the associated nonhomogeneous process at its maximum value, i.e., ${\widehat{\mathit{\lambda}}}_{i}={\mathit{\lambda}}_{i}\left(\mathrm{0}\right)$. We then generate the interarrival time to the next event t_{h} (${t}_{h}={t}^{*}{t}_{\mathrm{0}}$, where t^{*} is the time of occurrence of the event) from an exponential distribution with parameter $\mathit{\lambda}={\sum}_{n=\mathrm{1}}^{N}[{\mathit{\lambda}}_{n}{\mathbf{1}}_{\mathit{\{}n\notin \mathit{J}\mathit{\}}}+{\widehat{\mathit{\lambda}}}_{n}{\mathbf{1}}_{\mathit{\left\{}n\mathrm{\in}\mathit{J}\mathit{\right\}}}]$ (Step 6 in Fig. 5). The simulated event is then assigned either to one of the N hazard types ({H=h_{i}}) or to a null event from the nonhomogeneous Poisson process associated with the ith hazard type ($\mathit{\{}H={n}_{i}\mathit{\}},i\mathrm{\in}\mathit{J}$) (Step 7 in Fig. 5). The probability that the next event is the occurrence of the ith hazard type ({H=h_{i}}) is
and the probability that the next event to occur is a null event from the nonhomogeneous Poisson process associated with the ith hazard type ($\mathit{\{}H={n}_{i}\mathit{\}},i\mathrm{\in}\mathit{J}$) is
If the simulated event is a hazard event, Phases 1–3 described in Sect. 3.3.1 are performed (with the possible introduction of additional, nonhomogeneous Poisson processes). If the simulated event is a null event, it is discarded from the hazard event set. In both cases, we set ${t}_{\mathrm{0}}={t}^{*}$, and we update the rate of the nonhomogeneous Poisson process(es) to ${\widehat{\mathit{\lambda}}}_{i}={\mathit{\lambda}}_{i}({t}_{\mathrm{0}}{\widehat{t}}_{\mathrm{0}})$, where ${\widehat{t}}_{\mathrm{0}}$ is the time of occurrence of the event that started the nonhomogeneous Poisson process. Given the decaying nature of the rates of the nonhomogeneous processes considered herein, we assume that the effects of the altering interactions are forgotten whenever the difference between the updated rate of the process ${\widehat{\mathit{\lambda}}}_{i}$ and the rate of the process λ_{i} from the original rate curve falls below a prespecified threshold ε (i.e., if ${\widehat{\mathit{\lambda}}}_{i}{\mathit{\lambda}}_{i}<\mathit{\epsilon}$). The value of ε depends on the specific hazard, and it is subjective. Generally, ε can be selected based on experience or engineering judgment. Whenever ${\widehat{\mathit{\lambda}}}_{i}{\mathit{\lambda}}_{i}<\mathit{\epsilon}$ in the sequential MC procedure, we remove i from J, and we revert to using the original rate curve. A flowchart in Appendix A details the described sequential MC approach, with reference to each step shown in Fig. 5.
3.3.3 Incorporating slowonset events
This section details how to modify the procedure in Sect. 3.3.1 to incorporate possible slowonset events. Figure 6 visualizes the described procedure. The simulation of the ith slowonset event is herein modeled with a twostep approach. First, the event’s start is simulated using a homogeneous Poisson process with a distinct start rate ${\mathit{\lambda}}_{i}^{\mathrm{s}}$. Then, the start rate is replaced by an end rate ${\mathit{\lambda}}_{i}^{\mathrm{e}}$, which is used to model the conclusion of the event (with a corresponding homogeneous Poisson process). The rate change could be interpreted as an altering interaction of the hazard with itself. This method facilitates a realistic representation of events characterized by varying duration. The severity of such events is assumed to be nonvarying throughout the event. However, the procedure could be modified to include events with timevarying severity. Slowonset events whose presence could also affect the rate of additional processes (e.g., a drought affecting the rates of wildfires and floods) can be modeled by assigning a successive altering interactions to the start of the event (e.g., the start of the slowonset event causes a change in the rates associated with the secondary hazard types) and by associating the end of the event with the loss of memory of the interactions (i.e., the end of the slowonset event effectively acts as a memory loss event; see Sect. 3.3.1). As the memory loss is marked by the end of the slowonset event, no additional memoryloss event needs to be added in this case. To aid in the interpretation of Fig. 6, the reader may interpret hazard type 1 as droughts (slowonset), hazard type 2 as floods (rate decreases during droughts), and hazard type 3 as wildfires (rate increases during droughts).
A flowchart in Appendix A details the described sequential MC approach, with reference to each step shown in Fig. 6.
We now showcase the life cycle hazard event set simulations using the sequential MC method detailed in the previous sections. Figure 7 shows the hazard types considered in the example and their interaction. Figure 7 also reports the severity measure(s) adopted for the considered hazard types. There are no concurrent interactions considered in this example. However, the hazard type “heavy rain” is associated with multiple severity measures (i.e., intensity and duration), and the modeling of single hazards with multiple associated severity measures is equivalent to the modeling of concurrent hazards (see Sect. 3).
Earthquakes have been separated into mainshocks and aftershocks for modeling purposes. This distinction allows us to separate the rate curves for the two hazard types, with aftershocks having a rate equal to 0 before the occurrence of a mainshock and a conditional rate curve after the occurrence of a mainshock. The distinction also allows the system to retain a memory of the mainshock after the occurrence of the aftershocks. Without this distinction, the simulated occurrence of the first aftershock would redefine the rates of subsequent aftershocks, and the effects of the mainshock would be forgotten. More sophisticated models where the occurrence of early aftershocks in the sequence affects the rate of subsequent aftershocks can also be found in the literature (Ogata, 1998) but are not selected here for illustrative purposes. In the following, we report the reference for each of the adopted rate curves and surfaces for single hazards and each of the interaction models (in this order). Such references are also summarized in Table 1.
Iervolino et al. (2018)Tang and Cheung (2011)Yeo and Cornell (2009)Iervolino et al. (2018)Parker et al. (2015)Liu and Wang (2022)Samia et al. (2017)We note that the curves, surfaces, and/or models used in this numerical example have been developed for different locations worldwide. As such, the event sets obtained in this example shall not be associated with any specific location. The example showcases the potential of the proposed simulation method and shows that the required information can be retrieved from the literature. The collection of information for a specific location (which would require a tailored investigation and/or literature review) is outside the scope of the paper.
Mainshocks/aftershocks. The rate curve for the occurrence of the mainshock events, λ_{m}(m_{m}), is obtained from Iervolino et al. (2018). Namely, we use the exceedance curves for Zone 923 of the Italian earthquake hazard model, which corresponds to the L’Aquila region. The severity measure of earthquakes is expressed in terms of moment magnitude M_{w} (M_{m} for the mainshocks and M_{a} for the aftershocks). The minimum severity measure for both mainshocks and aftershocks is assumed as ${m}_{\mathrm{m},\mathrm{min}}={m}_{\mathrm{a},\mathrm{min}}=\mathrm{4.45}$, which is slightly higher than the one in the original reference (4.15). This is to reduce the number of small earthquakes in the simulated life cycle hazard event set to improve the clarity of this illustrative application. Aftershocks cannot occur independently from a mainshock. Therefore, their original rate curve is set to ${\mathit{\lambda}}_{\mathrm{a}}\left({m}_{\mathrm{a}}\right)=\mathrm{0}\forall {m}_{\mathrm{a}}$. All earthquakes are assumed to occur at the same location (or very close in space), such that the shortest Joyner–Boore distance (Joyner and Boore, 1981) of each earthquake from the location of interest is R_{jb}=20 km.
Heavy rain. The rate surface for the occurrence of heavy rain events, ${\mathit{\lambda}}_{r}({m}_{\mathrm{r},\mathrm{1}},{m}_{\mathrm{r},\mathrm{2}})$, is obtained from the intensity–duration–frequency curves reported by Tang and Cheung (2011) for the Hong Kong region in the GEO Report No. 261 (2011). The severity measures associated with the heavy rain events are the duration (M_{r,1}) expressed in hours (h) and the intensity (M_{r,2}) expressed in millimeters per hour (mm h^{−1}). The minimum considered severity measures are ${m}_{\mathrm{r},\mathrm{1},\mathrm{min}}=\mathrm{0.083}$ h (5 min) and ${m}_{\mathrm{r},\mathrm{2},\mathrm{min}}=\mathrm{0.893}$ mm h^{−1} for the duration and the intensity, respectively.
Landslide. We assume that landslide events cannot occur independently for this case study. They can only occur as secondary events of a successive interaction with earthquakes, heavy rain, and/or previous landslides. We also assume for simplicity that landslide events are not associated with any severity measure. The details of the abovementioned interactions are detailed below.
Mainshock–aftershock interactions. The interaction between mainshocks and aftershocks is modeled as a successive altering interaction following the modified Omori law (Yeo and Cornell, 2009). According to the Omori law, after the occurrence of a mainshock of magnitude m_{m} at time t_{m}, aftershocks occur following a nonhomogeneous Poisson process with a timevarying rate:
where t and t_{m} quantities are expressed in days, and a, b, c, and p are the parameters of the model. For this example, we assume (from Iervolino et al. (2018)) $a=\mathrm{1.66}$, b=0.96, c=0.03, and p=0.93. The severity of the aftershocks is simulated from the rate curve used for the mainshock (using Eq. 3).
Earthquake (mainshock/aftershock) and landslide interactions. The interaction among earthquakes and landslides is modeled as a successive triggering interaction following the model in Parker et al. (2015), calibrated based on data from the South Island region in New Zealand. The probability of a landslide given the occurrence of an earthquake with severity measure m_{w}, P_{0}(Lm_{w}), is
where SL is the local hillslope gradient (we assume SL =35° for the slope of this case study); NDS is the normalized distance from stream to ridge crest (we assume NDS =0.5 for this case study); c_{0}, c_{PGA}, c_{SL}, and c_{NDS} are parameters of the logistic regression used to fit Eq. (9) to the data; and PGA(m_{w},R_{jb}) is the peak ground acceleration caused by the earthquake at the location of interest. We obtain PGA(m_{w},R_{jb}) using the ground motion model in Huang and Galasso (2019), assuming that the soil type is rock and that the style of faulting of the earthquake is strikeslip fault:
where $\mathbf{b}=[{b}_{\mathrm{1}},{b}_{\mathrm{2}},{b}_{\mathrm{3}},{b}_{\mathrm{4}},{b}_{\mathrm{5}},{b}_{\mathrm{6}}]$ is a vector of unknown parameters. From the PGA model with spatial correlation in Huang and Galasso (2019), we select b_{1}=3.524, b_{2}=0.247, ${b}_{\mathrm{3}}=\mathrm{0.020}$, ${b}_{\mathrm{4}}=\mathrm{3.936}$, b_{5}=0.351, and b_{6}=12.417.
We note that the PGA(m_{w},R_{jb}) obtained from Eq. (10) is the median value at the location of interest, and a more refined analysis should consider the uncertainties associated with this value. However, a full probabilistic analysis of the ground motion is outside the scope of this example.
Heavyrain–landslide interactions. The interaction among heavy rain events and landslides is modeled as a successive triggering interaction following the model in Liu and Wang (2022). Namely, the probability of a landslide after a rainfall of intensity m_{r,2} is given by
where D_{c}(m_{r,2}) is the critical rainfall duration for slope instability associated with the intensity m_{r,2}. We use the critical rainfall duration for the slope after stabilization in Liu and Wang (2022), i.e.,
Landslide–landslide interactions. The interactions among subsequent landslides are modeled as a successive altering interaction following the model in Samia et al. (2017). The paper reports that the susceptibility of a slope to a landslide (s_{l}) increases by 15fold immediately after the occurrence of a previous landslide, and then it decreases exponentially over time with an exponential coefficient value of ${c}_{\mathrm{l}}=\mathrm{0.12}$. Hence, we define the timevarying susceptibility of the slope as
where t_{l1} is the time of occurrence of the first landslide event and ${t}_{{\mathrm{l}}^{*}}$ is the time of occurrence of the latest landslide event. Consequently, in the aftermath of a landslide event, we modify the probability of occurrence of landslides due to earthquakes and heavy rain events as $P\left(L\right{m}_{\mathrm{w}})={s}_{\mathrm{l}}(t\left){P}_{\mathrm{0}}\right(L\left{m}_{\mathrm{w}}\right)$ and $P\left(L\right{m}_{\mathrm{r},\mathrm{1}},{m}_{\mathrm{r},\mathrm{2}})={s}_{\mathrm{l}}(t\left){P}_{\mathrm{0}}\right(L{m}_{\mathrm{r},\mathrm{1}},{m}_{\mathrm{r},\mathrm{2}})$, respectively.
Figure 8 shows an example event set generated using the proposed method. The severity measure of the events is proportional to the diameter of the circles used to represent their occurrence. As expected, more severe mainshocks produce more severe aftershock sequences. The rate and severity measures (intensity and duration) of heavy rain events are generated based on their joint rate surface. Sequences of landslide events can also be found around years 17–25 in the system’s life cycle, triggered mostly by mainshock–aftershock sequences, and at year 40, triggered by heavy rain events.
Simulating a multihazard event set such as the one shown in Fig. 8 is computationally efficient since it amounts to sampling numbers from exponential distributions. Consequently, multiple realizations can be used to obtain relevant statistics for the system’s life cycle. For example, Table 2 reports the mean and median number of events, categorized by hazard type, during the system’s life cycle, obtained from 25 000 simulations. The difference between the mean and median value can serve as an effective indicator of the skewness of the associated distribution of events. When the mean and median are close, the distribution exhibits a scarcity of outliers, with realizations evenly distributed around the mean. Conversely, a notable difference between the mean and median signals heightened variability among realizations, suggesting the presence of multiple event sets with either a scarcity (e.g., zero) or a high number of events.
The simulated event sets may be used to identify the possible occurrence of independent hazard in close temporal proximity. To that purpose, Fig. 9 shows the mean (Fig. 9a) and median (Fig. 9b) number of pairwise hazard combinations throughout the life cycle of the system. A hazard combination is identified whenever two hazards occur with an interarrival time < 200 d (d signifies days), and no other hazard event occurs between them. The threshold of 200 d is arbitrary, and it has been selected solely for demonstration purposes. Nevertheless, it effectively represents a reasonable temporal interval that would pose challenges to the recovery from the initial hazard event before the occurrence of the subsequent hazard event (e.g., Opabola and Galasso, 2024). Heat maps such as the one in Fig. 9 can help identify hazard combinations that may pose significant challenges in terms of possible consequences on the physical assets of interest (which would require the modeling of Level II interactions). It can be observed that even hazard types that are not related by Level I interactions may occur close (in time) to each other. For example, there are, on average, 2.62 main shock events following the occurrence of a heavy rain event and 1.47 heavy rain events following a mainshock event, which might suggest that the interactions between such (independent) hazards might be of interest in a possible life cycle analysis of a structure placed in the investigated location. This sort of “coincidental” hazard combinations should not surprise the end users of the algorithm. In fact, such combinations have been observed on multiple occasions over the past century and are expected to increase due to climate change (e.g., Cutter et al., 2008). For example, typhoons were recorded in close temporal proximity to the great Kanto earthquake (Japan) of 1923 (e.g., Sasaki and Yamakawa, 2007) and the Hokkaido earthquake (Japan) of 2018 (e.g., Heidarzadeh et al., 2023). Although less frequent, these combinations are just as crucial as those influenced by causality (de Ruiter et al., 2020). An additional advantage of the proposed method is that it seamlessly integrates causal and coincidental event combinations within the same formulation.
Finally, joint distributions for the severity measures of multiple hazards can be obtained for each hazard combination of interest. For example, Fig. 10 shows the joint probability density function (obtained using the kernel density method) of the severity measures for the hazard combinations mainshock–mainshock (Fig. 10a), mainshock–aftershock (Fig. 10b), and mainshock–heavyrain (Fig. 10c). As expected, hazard combinations that are linked by Level I interactions show clear signs of correlation (mainshock–aftershock). In contrast, the joint probability density function (PDF) of the severity measures for independent hazards shows no apparent signs of correlation and can be approximated by the product of the marginal PDFs of the severity measures for the individual hazards (which can be obtained from Eq. 3).
The paper proposes a simple simulationbased approach to account for different types of hazard interactions in generating a life cycle multihazard event set, i.e., a sequence of events throughout a system’s life cycle and their associated characteristics and severities. We account for concurrent interactions, successive interactions where the secondary hazard event is immediately triggered by the primary one, and successive interactions where the primary hazard event affects the occurrence rate of the secondary one. Each interaction is incorporated into the simulation differently. Concurrent hazards are modeled based on the rate surface that defines the joint rate of the associated severity measures. Successive triggering interactions are incorporated through the conditional probability of occurrence of the secondary hazard type(s) and the conditional distribution of the associated severity measure. Successive altering interactions are modeled by modifying the secondary hazard type's rate curve. The different hazard events are modeled as a set of competing Poisson processes, which may be homogeneous or nonhomogeneous. The paper fills a gap in the literature for a quantitative interpretation of multihazard occurrences, translating the available, qualitative definitions and classifications into a systematic method to simulate event occurrences. The resulting simulation of one life cycle hazard event set is computationally efficient and can be repeated to obtain relevant statistics of hazard occurrences. By using competing Poisson processes and integrating within the same methodology both dependent and independent hazards, the proposed simulation method offers insights into not only the combination of hazards arising from causality (i.e., hazard interactions), but also those emerging from sheer temporal coincidence. The significance of these hazard combinations, especially in the context of their anticipated growth due to climate change effects in the coming years, should not be underestimated. By allowing the rate curves used as input to the model to be modified, the proposed algorithm allows us to incorporate such climate change aspects. The statistics of hazard occurrences can be used in analytical methods for life cycle analysis to obtain the associated statistics of impact/consequence metrics of interest for end users throughout the service life of a considered system. The simulated event sets can also be integrated into simulationbased frameworks for Level II interactions, i.e., impact/consequence interactions that can only occur through the exposed elements. This study specifically delves into the temporal dependence across hazards and does not explicitly discuss any spatial aspects/dependencies. The simulation of the hazard events in the example is based on the rate curves associated with their event characteristics (i.e., location and magnitude for earthquakes or rainfall intensity and duration for heavy rain events). A corresponding local intensity measure for earthquakes (i.e., the peak ground acceleration needed for landslide simulation) is obtained at a single location rather than across a region, although such an extension would be straightforward. In fact, local intensities (for each hazard event) are in general needed for a comprehensive analysis of Level II interactions. For analyses at the regional scale, in particular, the intensity measures are in the form of maps (hazard footprints) illustrating the hazard intensity measure's spatial variability across different locations, explicitly accounting for any spatial and crossintensity correlation (e.g., Jayaram and Baker, 2010a, b).
The code used to produce the results in this paper can be found in the following GitHub repository: https://github.com/LeandroIannacone/MultiHazardEventSetSimulation (last access: 2 May 2024) (https://doi.org/10.5281/zenodo.11061164, Iannacone, 2024).
No digital data repositories were accessed to produce the results in this work. Data from publicly available reports and papers are properly referenced in the text. In particular, Table 1 includes the reports and papers used for the case study.
LI: writing – original draft; conceptualization; methodology; formal analysis. KO: conceptualization; writing – review and editing. RG: supervision; conceptualization; writing – review and editing. CG: supervision; conceptualization; writing – review and editing.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
This article is part of the special issue “Methodological innovations for the analysis and management of compound risk and multirisk, including climaterelated and geophysical hazards (NHESS/ESD/ESSD/GC/HESS interjournal SI)”. It is not associated with a conference.
This research has been supported by UK Research and Innovation (grant no. NE/S009000/1) and by the European Union's Horizon Europe “Multihazard and riskinformed system for Enhanced local and regional Disaster risk management (MEDiate)” project under grant agreement no. 101074075 (UKRI reference number: 10049469 – Horizon Europe Guarantee).
This paper was edited by Avantika Gori and reviewed by two anonymous referees.
Abrahamson, N. A. and Bommer, J. J.: Probability and uncertainty in seismic hazard analysis, Earthq. Spectra, 21, 603–607, 2005. a
Adachi, T. and Ellingwood, B. R.: Serviceability of earthquakedamaged water systems: Effects of electrical power availability and power backup systems on system vulnerability, Reliab. Eng. Syst. Saf., 93, 78–88, 2008. a
Anagnos, T. and Kiremidjian, A. S.: A review of earthquake occurrence models for seismic hazard analysis, Probabilist. Eng. Mech., 3, 3–11, 1988. a
Apivatanagul, P., Davidson, R., Blanton, B., and Nozick, L.: Longterm regional hurricane hazard analysis for wind and storm surge, Coast. Eng., 58, 499–509, 2011. a
Cutter, S. L.: Compound, cascading, or complex disasters: what's in a name?, Environment: Science and Policy for Sustainable Development, 60, 16–25, 2018. a
Cutter, S. L., Barnes, L., Berry, M., Burton, C., Evans, E., Tate, E., and Webb, J.: A placebased model for understanding community resilience to natural disasters, Global Environ. Chang., 18, 598–606, 2008. a, b
Dehghani, N. L., Fereshtehnejad, E., and Shafieezadeh, A.: A Markovian approach to infrastructure lifecycle analysis: Modeling the interplay of hazard effects and recovery, Earthquake Eng. Struct. D., 50, 736–755, 2021. a, b
Der Kiureghian, A. and Ang, A. H.: A faultrupture model for seismic risk analysis, B. Seismol. Soc. Am., 67, 1173–1194, 1977. a
de Ruiter, M. C., Couasnon, A., van den Homberg, M. J., Daniell, J. E., Gill, J. C., and Ward, P. J.: Why we can no longer ignore consecutive disasters, Earth's Future, 8, e2019EF001425, https://doi.org/10.1029/2019EF001425, 2020. a, b
Di Baldassarre, G., Nohrstedt, D., Mård, J., Burchardt, S., Albin, C., Bondesson, S., Breinl, K., Deegan, F. M., Fuentes, D., Lopez, M. G., Granberg, M., Nyberg, L., Nyman, M. R., Rhodes, E., Troll, V., Young, S., Walch, C., and Parker, C. F.: An integrative research framework to unravel the interplay of natural hazards and vulnerabilities, Earth's Future, 6, 305–310, 2018. a
Douglas, J. and Edwards, B.: Recent and future developments in earthquake ground motion estimation, EarthSci. Rev., 160, 203–219, 2016. a
Fadhel, S., RicoRamirez, M. A., and Han, D.: Uncertainty of intensity–duration–frequency (IDF) curves due to varied climate baseline periods, J. Hydrol., 547, 600–612, 2017. a
Gardoni, P. and LaFave, J. M.: Multihazard approaches to civil infrastructure engineering: Mitigating risks and promoting resilence, Springer, https://doi.org/10.1007/9783319297132_1, 2016. a
Gentile, R., Cremen, G., Galasso, C., Jenkins, L. T., Manandhar, V., Menteşe, E. Y., Guragain, R., and McCloskey, J.: Scoring, selecting, and developing physical impact models for multihazard risk assessment, Int. J. Disast. Risk Re., 82, 103365, https://doi.org/10.1016/j.ijdrr.2022.103365, 2022. a
Gill, J. C. and Malamud, B. D.: Reviewing and visualizing the interactions of natural hazards, Rev. Geophys., 52, 680–722, 2014. a, b, c
Gill, J. C. and Malamud, B. D.: Anthropogenic processes, natural hazards, and interactions in a multihazard framework, EarthSci. Rev., 166, 246–269, 2017. a, b
Heidarzadeh, M., Miyazaki, H., Ishibe, T., Takagi, H., and Sabeti, R.: Field surveys of September 2018 landslidegenerated waves in the Apporo dam reservoir, Japan: combined hazard from the concurrent occurrences of a typhoon and an earthquake, Landslides, 20, 143–156, 2023. a
Huang, C. and Galasso, C.: Groundmotion intensity measure correlations observed in Italian strongmotion records, Earthqu. Eng. Struct. D., 48, 1634–1660, 2019. a, b
Iacoletti, S., Cremen, G., and Galasso, C.: Integrating Long and ShortTerm Time Dependencies in SimulationBased Seismic Hazard Assessments, Earth Space Sci., 9, e2022EA002253, https://doi.org/10.1029/2022EA002253, 2022. a, b
Iannacone, L.: Multi Hazard Event Set Simulation, Zenodo [code], https://doi.org/10.5281/zenodo.11061164, 2024.
Iervolino, I. and Giorgio, M.: Comment on “How well does Poissonian probabilistic seismic hazard assessment (PSHA) approximate the simulated hazard of epidemictype earthquake sequences?” by Shaoqing Wang, Maximilian J. Werner, and Ruifang Yu, B. Seismol. Soc. Am., 112, 2758–2761, 2022. a
Iervolino, I., Giorgio, M., and Polidoro, B.: Sequencebased probabilistic seismic hazard analysis, B. Seismol. Soc. Am., 104, 1006–1012, 2014. a
Iervolino, I., Chioccarelli, E., and Giorgio, M.: Aftershocks’ effect on structural design actions in Italy, B. Seismol. Soc. Am., 108, 2209–2220, 2018. a, b, c, d
Jayaram, N. and Baker, J. W.: Considering spatial correlation in mixedeffects regression and the impact on groundmotion models, B. Seismol. Soc. Am., 100, 3295–3303, 2010a. a, b
Jayaram, N. and Baker, J. W.: Efficient sampling and data reduction techniques for probabilistic seismic lifeline risk assessment, Earthqu. Eng. Struct. D., 39, 1109–1131, 2010b. a, b
Joyner, W. B. and Boore, D. M.: Peak horizontal acceleration and velocity from strongmotion records including records from the 1979 Imperial Valley, California, earthquake, B. Seismol. Soc. Am., 71, 2011–2038, 1981. a
Kappes, M. S., Keiler, M., von Elverfeldt, K., and Glade, T.: Challenges of analyzing multihazard risk: a review, Nat. Hazards, 64, 1925–1958, 2012. a
Lewis, P. W. and Shedler, G. S.: Simulation of nonhomogeneous Poisson processes by thinning, Nav. Res. Logist. Q., 26, 403–413, 1979. a, b, c
Liu, X. and Wang, Y.: Quantifying annual occurrence probability of rainfallinduced landslide at a specific slope, Comput. Geotech., 149, 104877, https://doi.org/10.1016/j.compgeo.2022.104877, 2022. a, b, c
Marzocchi, W., GarciaAristizabal, A., Gasparini, P., Mastellone, M. L., and Di Ruocco, A.: Basic principles of multirisk assessment: a case study in Italy, Nat. Hazards, 62, 551–573, 2012. a
Mignan, A., Wiemer, S., and Giardini, D.: The quantification of lowprobability–highconsequences events: part I. A generic multirisk approach, Nat. Hazards, 73, 1999–2022, 2014. a, b
Mignot, E. and Dewals, B.: Hydraulic modelling of inland urban flooding: recent advances, J. Hydrol., 609, 127763, https://doi.org/10.1016/j.jhydrol.2022.127763, 2022. a
Neri, A., Aspinall, W. P., Cioni, R., Bertagnini, A., Baxter, P. J., Zuccaro, G., Andronico, D., Barsotti, S., Cole, P. D., Ongaro, T. E., Hinks, T. K., Macedonio, G., Papale, P., Rosi, M., Santacroce, R., and Woo, G.: Developing an event tree for probabilistic hazard and risk assessment at Vesuvius, J. Volcanol. Geoth. Res., 178, 397–415, 2008. a, b
Nofal, O. M., Amini, K., Padgett, J. E., van de Lindt, J. W., Rosenheim, N., Darestani, Y. M., Enderami, A., Sutley, E. J., Hamideh, S., and DuenasOsorio, L.: Multihazard sociophysical resilience assessment of hurricaneinduced hazards on coastal communities, Resilient Cities and Structures, 2, 67–81, 2023. a
Ogata, Y.: Spacetime pointprocess models for earthquake occurrences, Ann. I. Stat. Math., 50, 379–402, 1998. a
Opabola, E. A. and Galasso, C.: A probabilistic framework for postdisaster recovery modeling of buildings and electric power networks in developing countries, Reliab. Eng. Syst. Safe., 242, 109679, https://doi.org/10.1016/j.ress.2023.109679, 2024. a
Otárola, K., Iannacone, L., Gentile, R., and Galasso, C.: A Markovian framework for multihazard lifecycle consequence analysis of deteriorating structural systems, in: Proceedings of the 14th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP14, Dublin, Ireland, 9–13 July 2023, https://discovery.ucl.ac.uk/id/eprint/10180497 (last access: 2 May 2024), 2023. a, b
Otárola, K., Iannacone, L., Gentile, R., and Galasso, C.: Multihazard lifecycle consequence analysis of deteriorating engineering systems, Struct. Saf., in review, 2024. a, b
Pagani, M., GarciaPelaez, J., Gee, R., Johnson, K., Poggi, V., Silva, V., Simionato, M., Styron, R., Viganò, D., Danciu, L., Monelli, D., and Weatherill, G.: The 2018 version of the global earthquake model: Hazard component, Earthqu. Spectra, 36, 226–251, 2020. a
Parker, R. N., Hancox, G. T., Petley, D. N., Massey, C. I., Densmore, A. L., and Rosser, N. J.: Spatial distributions of earthquakeinduced landslides and hillslope preconditioning in the northwest South Island, New Zealand, Earth Surf. Dynam., 3, 501–525, https://doi.org/10.5194/esurf35012015, 2015. a, b, c
Pescaroli, G. and Alexander, D.: Understanding compound, interconnected, interacting, and cascading risks: a holistic framework, Risk Anal., 38, 2245–2257, 2018. a
Ritschel, C., Ulbrich, U., Névir, P., and Rust, H. W.: Precipitation extremes on multiple timescales – Bartlett–Lewis rectangular pulse model and intensity–duration–frequency curves, Hydrol. Earth Syst. Sci., 21, 6501–6517, https://doi.org/10.5194/hess2165012017, 2017. a
Sadegh, M., Moftakhari, H., Gupta, H. V., Ragno, E., Mazdiyasni, O., Sanders, B., Matthew, R., and AghaKouchak, A.: Multihazard scenarios for analysis of compound extreme events, Geophys. Res. Lett., 45, 5470–5480, 2018. a
Samia, J., Temme, A., Bregt, A., Wallinga, J., Guzzetti, F., Ardizzone, F., and Rossi, M.: Characterization and quantification of path dependency in landslide susceptibility, Geomorphology, 292, 16–24, 2017. a, b
Sasaki, H. and Yamakawa, S.: Natural hazards in Japan, in: International perspectives on natural disasters: occurrence, mitigation, and consequences, Springer, 163–180, https://doi.org/10.1007/9781402028519_8, 2007. a
Selva, J.: Longterm multirisk assessment: statistical treatment of interaction among risks, Nat. Hazards, 67, 701–722, 2013. a, b, c, d
Shaluf, I. M.: An overview on disasters, Disaster Prev. Manag., 16, 687–703, 2007. a
Tang, C. and Cheung, S.: Frequency analysis of extreme rainfall values, in: Proceedings of the HKIE Geotechnical Division Annual Seminar on Landslide Risk Reduction through Works: ThirtyFive Years of Landslip Preventive Measures Programme and Beyond, Hong Kong Institution of Engineers Hong Kong, https://www.cedd.gov.hk/eng/publications/geo/georeports/geo_rpt261/index.html (last access: 2 May 2024), 2011. a, b
UNISDR: (United Nations International Strategy for Disaster Reduction) Hyogo framework for action 2005–2015: Building the resilience of nations and communities to disasters, Geneva: UNISDR, https://www.undrr.org/publication/hyogoframeworkaction20052015buildingresiliencenationsandcommunitiesdisasters (last access: 2 May 2024), 2005. a, b
Utsu, T.: Aftershocks and earthquake statistics (1): Some parameters which characterize an aftershock sequence and their interrelations, Journal of the Faculty of Science, Hokkaido University, Series 7, Geophysics, 3, 129–195, 1970. a, b
Wang, Y. and Rosowsky, D. V.: Characterization of joint wind–snow hazard for performancebased design, Struct. Saf., 43, 21–27, 2013. a
Westcott, M.: A note on the nonhomogeneous Poisson cluster process, J. Appl. Probab., 14, 396–398, 1977. a
Yeo, G. L. and Cornell, C. A.: A probabilistic framework for quantification of aftershock groundmotion hazard in California: Methodology and parametric study, Earthqu. Eng. Struct. D., 38, 45–60, 2009. a, b
Zaghi, A. E., Padgett, J. E., Bruneau, M., Barbato, M., Li, Y., MitraniReiser, J., and McBride, A.: Establishing common nomenclature, characterizing the problem, and identifying future opportunities in multihazard design, J. Struct. Eng., 142, H2516001, https://doi.org/10.1061/(ASCE)ST.1943541X.0001586, 2016. a, b, c, d, e
Zscheischler, J., Westra, S., Van Den Hurk, B. J., Seneviratne, S. I., Ward, P. J., Pitman, A., AghaKouchak, A., Bresch, D. N., Leonard, M., Wahl, T., and Zhang, X.: Future climate risk from compound events, Nat. Clim. Change, 8, 469–477, 2018. a