Simulating Multihazard Interactions Using Higher-Order Network Analysis

Multihazard risk analysis requires appropriate consideration of two levels of interactions—Level I (or direct) and Level II (or indirect) interactions (Zaghi et al. 2016). Level I interactions occur due to the inherent interrelationships of hazards and are independent of the presence of physical components. Level II interactions look at hazard interrelationships through site effects and impacts on physical components, systems, and networks (e.g., flood occurrence due to the earthquake-induced collapse of a dam or pipelines, wildfire occurrence from earthquake-induced damage to gas refineries, or landslide occurrence due to the earthquake-induced collapse of a retaining structure). The framework presented in this study focuses only on Level I interactions. Future studies will explore extending the methodology to consider Level II interactions.

Appropriate multihazard analysis of Level I interactions entails identifying all potential hazard sources in a geographical region of interest, defining and characterizing their potential interactions, defining relevant occurrence models for each hazard considering such interactions, and probabilistically simulating multihazard scenarios over a defined space-time domain. In this study, hazard sources collectively refer to various phenomena (natural or human-induced) that result in events that can potentially lead to disasters. Sometimes, these hazard sources are physically tied to specific locations. For example, earthquakes are linked to geological faults, and riverine flooding is associated with rivers. A flowchart for multihazard analysis capturing the steps mentioned previously is presented in Fig. 1. The subsequent subsections provide a detailed description of each step.

Hazard analysis characterizes the occurrence time, location, severity, spatial and temporal scales, and (local) intensities of potential hazards in a region of interest. A fundamental understanding of relevant hazard sources is critical in multihazard analysis. This entails, for instance, collecting information on the location and rupture history of seismic faults, volcanoes, information relevant to flood phenomena (e.g., rainfall and other meteorological data, river/stream locations, and streamflow and water level data), hurricane phenomena, and liquefaction- and landslide-prone areas (e.g., soil characteristics). Also, details on tsunamigenic processes in the region of interest must be collated in tsunami-prone areas.

Once all the hazard sources have been identified, all potential hazard interactions in the region must be defined. A wide range of terminologies are adopted in multihazard risk assessment to classify hazard interactions [see Kappes et al. (2010), Gill and Malamud (2014), Liu et al. (2016), Zaghi et al. (2016), and Tilloy et al. (2019)]. We note that it is outside the scope of this study to validate (or otherwise) any of the published terminologies. However, capturing hazard interactions defined in the existing literature is essential for a holistic multihazard modeling/simulation framework, as proposed in this study.

A review of the definitions discussed in the aforementioned literature suggests that the hazard interaction classifications discussed in Tilloy et al. (2019) seem to encompass well the classifications presented by Kappes et al. (2010), Gill and Malamud (2014), and Liu et al. (2016). The Zaghi et al. (2016) classifications for interacting hazards (i.e., successive and concurrent) focus on the time interval between a primary and a secondary hazard.

The framework presented in this paper has been developed to capture the definitions presented by Tilloy et al. (2019) and Zaghi et al. (2016). The independent (or coincidental) hazard interaction helps capture the spatial and temporal overlapping of the impact(s) of two or more independent hazards (i.e., no triggering relationship). An example is the joint occurrence of wind and snow hazards. The triggering hazard interaction captures a case where a primary hazard triggers a secondary hazard (e.g., earthquake-tsunami or earthquake-liquefaction sequence) from the primary hazard. An alteration hazard interaction captures a case where the occurrence of a primary hazard can alter the dispositions (i.e., occurrence rate and severity) of a secondary hazard. An example of an alteration interaction is a wildfire destroying vegetation and increasing the local severity of any flood occurrence. The compound interaction characterizes a case where the co-occurrence of two different hazards can result from the same large-scale process. An example is the co-occurrence of sea surge and river flooding (both from a tropical cyclone) to cause compound flooding. The Tilloy et al. (2019) definition of the mutual exclusion interaction encompasses negative dependence (e.g., flood and fire), which can be mathematically expressed as altered severity and occurrence rate (i.e., alteration interaction).

It is also important to identify the hazards based on their temporal scale. Depending on the temporal scale over which they occur, hazards can be categorized into rapid-onset (also called sudden-onset) and slow-onset processes (e.g., Jones 1993). Rapid-onset processes are sudden and unexpected events that have an immediate impact. Their duration is typically short, ranging from tens of seconds, as in the case of earthquakes, to several hours, as seen with volcanic eruptions and flash floods. Conversely, slow-onset processes extend over a longer period, from days to weeks for short-term events like hurricanes, river floods, and wildfires, to months or even years for long-term events such as droughts, famines, and pandemics (Jones 1993).

This study adopts directed hypergraphs in modeling the complex interactions of hazards. Hypergraphs provide the most unconstrained description of higher-order interactions (e.g., Berge 1989). A hypergraph—$H(V,E)$—is a generalization of a graph in which an edge can join any number of vertices. In contrast, an edge connects exactly two vertices in an ordinary graph. Specifically, a hypergraph is defined by a set of vertices $V=\{v_i|i=1,\dots ,p\}$ joined by a set of edges $E=\{E_j|j=1,\dots ,q\}$ that specify which vertices participate in an interaction. Hypergraph theory has been applied in various fields, e.g., chemistry (e.g., Konstantinova and Skorobogatov 2001), image processing (e.g., Bretto and Gillibert 2005), and social media (e.g., Amini et al. 2022), to study complex networks.

The proposed framework uses the hypergraph theory to describe a multihazard-prone region of hazard vertices and edges. The hazard vertices represent the possible hazard phenomena—e.g., earthquake from faults, river flooding, hurricane, and drought. Each hazard vertex has its attributes (e.g., occurrence rate, meteorological parameters, geotechnical parameters, or rupture characteristics). Hazard interactions are represented as hazard edges in the hypergraph. Each edge has one primary and one or more secondary hazards. Herein, we define the primary hazard for an edge as a hazard that can occur without being triggered. Secondary hazards in an edge can have their occurrence triggered or altered by the primary hazard in the edge. To further explain this definition, Fig. 2(a) shows an example of a map of hazard sources in a hypothetical region with two seismic faults (Fault 1 and Fault 2), a river, a liquefiable site, and an ocean with tsunamigenic zones. For the region shown in Fig. 2(a), let us assume mainshock earthquake events from Fault 1 ($h_1$) can trigger aftershocks in the same fault and also trigger liquefaction in the liquefiable soil ($h_4$)—i.e., Edge $E_1$ [see hypergraph in Fig. 2(b)]. Edge $E_2$ is characterized by mainshock earthquake events from Fault 2 ($h_2$) that can trigger aftershocks, liquefaction in the liquefiable soil ($h_4$) and tsunamis in the ocean ($h_5$). Also, let us assume the tsunami triggered by an earthquake in Fault 2 can travel up the river and low-lying areas along the riverbank to cause flooding ($h_3$). We further assume that background meteorological and climatic conditions (e.g., rainfall) can cause fluvial flooding ($h_3$) in the region—Edge $E_3$. Fig. 2(b) shows the five hazard vertices in the hypothetical region in Fig. 2(a), representing the interactions in the hypothetical region. As shown in Fig. 2(b), each primary hazard and the links between the hazards in each edge are represented by occurrence models. Further discussions on occurrence models are presented subsequently in this paper.

For directed hypergraphs, in each edge, the vertex can further be divided into source (so-called in hypergraph theory and not related to hazard sources as used in this paper) vertex set and target vertex set. For our purpose, the source vertex set represents the primary hazards in each edge, and the target vertex set represents all secondary hazards linked to each primary hazard for each edge. Eqs. (1) and (2) show the source ($H_s$) and target ($H_t$) matrices for Fig. 2(b). The number of edges ($j=1,\dots ,q$) in Eqs. (1) and (2) corresponds to the number of primary hazards. The number of rows ($i=1,\dots ,p$) in the source ($H_s$) and target ($H_t$) matrices corresponds to the total number of hazard sources in the region. For computation ease, the primary hazards are indexed before the secondary hazards, such that an identity matrix for the primary hazards can be derived from $H_s$ [see Eq. (1)]—i.e., $H_s(i,j)=1$, if $i=j$ or zero otherwise. Also, $H_t$ $(i,j)=1$ if the secondary hazard $h_i$ has a relationship/dependence with the primary hazard $h_j$ (in hypergraph $E_j$) or zero otherwise. In cases where a primary hazard triggers or alters the occurrence of a secondary hazard from the same source (e.g., mainshock-aftershock sequence from the same fault), the hazard is represented in both $H_s$ and $H_t$—i.e., $H_s(i,j)=H_t(i,j)=1$. We can also represent the primary hazards as a source unit column vector $H_{s,\mathrm{col}}$, with size ($q\times 1$)where $h_1$ = earthquake from Fault 1; $h_2$ = earthquake from Fault 2; $h_3$ = flooding from river; $h_4$ = liquefaction from liquefiable site; and $h_5$ = tsunami from ocean.

The incidence matrix $\widehat{H}$ with size ($p\times q$) for a directed hypergraph is given as follows:where $T(E)$ and $Q(E)$ = target and source vertices for edge $E$, respectively.

The directed hypergraph in Fig. 2(b) can be represented as shown in Eq. (4). In cases where a primary hazard triggers or alters the occurrence of a secondary hazard from the same source, the corresponding element in the incidence matrix is equal to unity—i.e., $\widehat{H}(i,j)=1$

The hypergraph shown in Fig. 2(b) is static (i.e., with no temporal attributes). Hence, the incidence matrix $\widehat{H}(V,E)$ developed from the static hypergraph [see Eq. (4)] cannot characterize the stochastic nature of multihazard interactions over a given time range. Such complex interactions and their stochastic nature cannot be represented with a static hypergraph but a temporal hypergraph. A temporal hypergraph $\widehat{H}(V,E,T)$ is an assemble of hypergraphs $\widehat{H}(V,E)$ with temporal attributes (e.g., time of hazard occurrence) associated with the vertices or edges—i.e., the temporal hypergraph $\widehat{H}(V,E,T)$ is an assembly of $m$ static hypergraphs ${\widehat{H}}_l(V,E)$ over a period [$0,T$] (where $l=1,\dots ,m$). The temporal attributes of the hypergraph can be defined using occurrence models and probabilistic simulation of events—described later in this section.

Reliable estimation of the frequency, timing, severity, and duration (especially for slow-onset events) of hazards is important for multihazard scenario simulation. The term severity is used in this study to refer to both event characteristics and local intensity measures. In fact, some aspects of hazard interactions 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 intensity at the slope location).

A set of occurrence and severity (including severity-time models for slow-onset hazards) models to capture the occurrence and interactions of hazards in the hypergraph edges must be defined for multihazard scenario simulation. The choice of occurrence and severity models depends on factors such as the nature of the hazard, available data, modeling objectives, and computational resources. Providing a detailed overview of all occurrence and severity models is outside the scope of this study. However, we highlight a few occurrence and severity models to support the discussions in this paper.

Hazard occurrence can be represented using homogeneous (HPP) and nonhomogeneous (NHPP) Poisson processes [e.g., hurricane—(Elsner et al. 2004); earthquakes—(Kagan and Jackson 1991); wind—(Gusclla 1991); floods—(Shane and Lynn 1964)], Weibull processes [e.g., flood—(Heo et al. 2001); earthquakes—(Hagiwara 1974); wind—(Xiao et al. 2006)], Brownian passage time (BPT) models [e.g., earthquakes—(Matthews et al. 2002)], among many others. Occurrence models typically combine various models that are based on input features such as environmental variables, land use, and other socioeconomic parameters. For example, flood occurrence can be simulated by combining rainfall-run-off models and hydraulic models (Rasouli et al. 2012; Schumann et al. 2013; Teng et al. 2017; Gori et al. 2019; Moftakhari et al. 2019). Studies (e.g., Gabriel et al. 2017; Oliveira et al. 2012) have also developed wildfire occurrence models that link meteorological factors (e.g., temperature, precipitation, wind speed) and land use. Also, drought occurrence can be simulated using the standardized precipitation index (SPI) (McKee et al. 1993) and/or the standardized soil moisture index SSMI (Hao and AghaKouchak 2013).

For triggering and altering interactions, the occurrence probability (or the rate curve relating the event characteristics to their corresponding exceedance rates) of a secondary hazard (say $h_{\sec }$) over a time interval ($0,t_{\mathrm{int}}$] is typically conditioned on the severity ($S_{h,\mathrm{pri}}$) and the occurrence rate (${\lambda }_{h,\mathrm{pri}}$) of the primary hazard (say $h_{\mathrm{pri}}$). For example, for an HPP, assuming that $h_{\sec }$ is triggered if $S_{h,\mathrm{pri}}$ reaches a certain intensity threshold $s^{*}$, the probability of occurrence of $h_{\sec }$ is given by Eq. (5). For example, as demonstrated in the case study presented subsequently in this paper, the simulation of the occurrence of earthquake-induced liquefaction could entail combining probabilistic liquefaction index models (e.g., Idriss and Boulanger 2010) with probabilistic local seismic intensity measure models [e.g., Boore et al. (2014) model] and earthquake occurrence models (e.g., Cornell and Vanmarcke 1969)

The alteration interaction (i.e., the increased or decreased occurrence rate or severity of a secondary hazard) can generally be expressed as a function of the occurrence rate ${\lambda }_{h,\mathrm{pri}}$, severity of the primary hazard $S_{h,\mathrm{pri}}$, the conditional occurrence model of the secondary hazard, ${\lambda }_{h,\sec }^o(t)$—e.g., the altered occurrence rate of a secondary hazard can be expressed as ${\lambda }_{h,\sec }(t)=f({\lambda }_{h,\sec }^o(t),S_{h,\mathrm{pri}},{\lambda }_{h,\mathrm{pri}})$. Alteration interactions can be captured using a single or a combination of models. Several studies (e.g., Liu et al. 2010; Peters and Iverson 2017; Taufik et al. 2017) have developed models to simulate wildfire occurrence amplification in drought-affected areas. Also, climate change models (e.g., Lee et al. 2021) can be linked to occurrence models for hydrometeorological hazards to capture altering relationships between climate change and the meteorological factors present in these occurrence models. There may be instances where the relationship between two hazards is not straightforward (e.g., studies have considered mainshock-aftershock sequences as either triggered or altered relationships). In such situations, risk analysts can choose the occurrence models that best align with the relationship they assume.

For rapid-onset processes, only single-point interarrival times, simulated from the defined stochastic process, are needed for scenario modeling. For example, for an NHPP, a thinning algorithm (Lewis and Shedler 1979) can be used to generate the interarrival times and the event duration is assumed to be zero (i.e., negligible with respect to the typical time range—e.g., 50 or 100 years) for multihazard scenario simulation. Similarly, severity models for rapid-onset processes can be used to generate a single point severity metric [e.g., tsunami inundation depth—see Abe (1995) and Kulikov et al. (2005)]. For slow-onset processes, the process duration and severity fluctuation during the process need to be characterized. For example, studies (e.g., Kirono et al. 2020; Sharma and Panu 2012) provide probabilistic models for predicting the duration of drought events as a function of SPI. The severity models for slow-onset typically generate a severity-time series—e.g., time series of drought index (e.g., Kirono et al. 2020).

Once the occurrence and severity models are defined for each hazard type of interest, a probabilistic simulation of multihazard events can be carried out. The flowchart for the multihazard scenario simulation within a desired time interval [$0,T$] is presented in Figs. 3 and 4. To start the analysis, we initialize the initial time $t_0=0$ and set ${\widehat{H}}_0=\mathrm{Ø}$, and $S_0=\mathrm{Ø}$. Subsequently, the arrival times $t_{hp}^{*}(j)$ of all primary hazards in the source column vector ($H_{s,\mathrm{col}}$) are simulated using the defined occurrence models. We then set the occurrence times $[t_{\mathrm{occ},hp}(j)]$ as the sum of $t_{hp}^{*}(j)$ and $t_o$. Subsequently, we set the current time $t_{\mathrm{occ},\text{curr},1}$ to be the minimum of the simulated occurrence times of the primary hazards [i.e., $t_{\mathrm{occ},\text{curr},1}=t_{\mathrm{occ},hp}(\underset{j}{\mathrm{argmin}}{t}_{\mathrm{occ},hp}(j))$]. As shown in Fig. 3, if $t_{\mathrm{occ},\text{curr},1}$ exceeds the defined $T$, we assume that no event will happen within the defined time interval [$0,T$] and end the simulation. Otherwise, we proceed with simulating the occurrence of the current primary hazard $h_{p,\mathrm{curr}}$ at the current time. The hypergraph incidence matrix at time $t_{\mathrm{occ},\text{curr},1}$ (i.e., ${\widehat{H}}_{t_{\mathrm{occ},\text{curr},1}}$) is updated to reflect the occurrence of $h_{p,\mathrm{curr}}$ (i.e., the hypergraph at the time of occurrence of the primary hazard is updated accordingly—${\widehat{h}}_{i,j}$ equals 1 for simulated primary hazard). We then use the defined severity models to simulate the severity of $h_{p,\mathrm{curr}}$ at $t_{\mathrm{occ},\text{curr},1}$, $S_{t_{\mathrm{occ},\text{curr},1}}$. Fig. 5 illustrates a hypothetical multihazard analysis for the region in Fig. 2 to generate an array of incidence and severity matrices across the time interval.

The next step is to simulate the temporal features of $h_{p,\mathrm{curr}}$—i.e., the event duration $t_{\mathrm{dur},\mathrm{curr}}$ and severity-time series. For rapid-onset processes (e.g., earthquakes and tsunami), event duration can be assumed to be zero, and the event is represented as a point in time—i.e., event end time $t_{\mathrm{occ},\text{curr},2}=t_{\mathrm{occ},\text{curr},1}$ and only appear in one hypergraph incidence matrix and severity matrix. For example, Fig. 5 depicts $h_1$, $h_2$, $h_4$, and $h_5$ as rapid-onset hazards; hence, each of their occurrences is represented in a single incidence matrix and a single severity matrix.

For slow-onset processes, $t_{\mathrm{occ},\text{curr},2}=t_{\mathrm{occ},\text{curr},1}+t_{\mathrm{dur},\text{curr}}\ge t_{\mathrm{occ},\text{curr},1}$. The associated severity fluctuation between times $t_{\mathrm{occ},\mathrm{curr},1}$ and $t_{\mathrm{occ},\mathrm{curr},2}$ for slow-onset processes can be simulated using relevant duration forecasting approaches. Slow-onset processes can appear in more than one hypergraph and severity incidence matrices. The number of defined hypergraph and severity incidence matrices depends on a predefined timestep. The timestep should be set such that the peaks and troughs of the severity-time series are properly captured. This is crucial for properly defining the corresponding severity matrices within the time range (i.e., $t_{\mathrm{occ},\mathrm{curr},1}$ and $t_{\mathrm{occ},\mathrm{curr},2}$). Fig. 5 depicts $h_3$ as a slow-onset process occurring between time $t=b$ and $t=e$ with a simulated severity-time history (e.g., maximum inundation depth versus time). The corresponding severity ($s_{i,j}^t$—where, $i$ = hazard, $j$ = edge, $t$ = current time) at each timestep can be derived from the severity-time series.

The next step is to simulate the impact of the occurrence of the current primary hazard on other hazards. As shown in Fig. 3, if the current primary hazard has no interactions with other hazards, we proceed to simulate a new interarrival time $t_{\mathrm{int},\mathrm{curr}}^{*}$ and next occurrence time of the primary hazard $[t_{\mathrm{occ},hp}(\mathrm{curr})]$, and then update $t_{\mathrm{occ},hp}$. Subsequently, the entire hazard occurrence simulation process is repeated with the updated for the primary hazard with the updated $t_{\mathrm{occ},hp}$.

As shown in Fig. 3, if the current primary hazard has an altered or triggered interaction with other hazards, we simulate each interaction using the flowchart presented in Fig. 4. If an altered interaction exists between the simulated primary hazard and other hazards, for each altered hazard, we update their corresponding occurrence and severity models accordingly (i.e., using the predefined occurrence and severity models). If the altered hazard is contained in $H_{s,\mathrm{col}}$, we use the updated occurrence model to simulate a new arrival time $t_{hp}^{*}(\mathrm{alt})$ and new occurrence time $t_{\mathrm{occ},hp}(\mathrm{alt})=t_{hp}^{*}(\mathrm{alt})+t_{\mathrm{occ},\mathrm{curr},1}$ for this altered primary hazard and update $t_{\mathrm{occ},hp}$.

In the case of triggered interactions, we simulate the occurrence time and severity of the triggered hazards (including temporal features) using corresponding conditional occurrence and severity models. Similar to the process for the primary hazard, the start time $t_{\mathrm{occ},\mathrm{trig},\mathrm{curr},1}$ and end time $t_{\mathrm{occ},\mathrm{trig},\mathrm{curr},2}$ (if slow-onset process) of the triggered hazard is simulated. Furthermore, we update the hypergraph incidence matrices and the severity matrices–i.e., ${\widehat{h}}_{i,j}$ for ${\widehat{H}}_t$ equals $-1$ for the simulated secondary hazard (as a result of the occurrence of the primary hazard) at time $t$ (see Fig. 5). As shown in Fig. 4, the framework also captures scenarios where a triggered hazard could alter the severity and/or occurrence models of other hazards.

As shown in Fig. 3, the simulation continues until a $t_{\mathrm{occ},\text{curr},1}\le T$ is not achieved. The simulation of multihazard events will generate a number of hypergraph incidence matrices similar to those in Fig. 5. Once the simulation is complete [i.e., $t_{\mathrm{occ},hp}(j)>T$], all the hypergraphs can be aggregated together to develop the multihazard scenario shown in Fig. 5.

The identification of concurrent independent hazards is important for multihazard risk management purposes. For example, about 70 countries reported flooding events during the COVID-19 pandemic (Simonovic et al. 2021). Apart from overlapping independent occurrences, there is also a coincidence window for the joint impact of independent hazards occurring in close proximity. For example, a catastrophic flood occurred roughly a month after the 2023 Turkish-Syria earthquake sequence, further affecting the earthquake victims (Kirby 2023).

The joint occurrence of independent hazards can be captured using event coincidence analysis (ECA) (e.g., Donges et al. 2016). This study adopts the sliding window procedure for temporal hypergraphs, widely used for signal processing (Myers et al. 2023), for event coincidence analysis (i.e., to capture the concurrent occurrence of independent hazards). In the sliding window procedure, a window of specified width (hereafter referred to as the coincidence interval) $w$, slides along the time domain of a temporal hypergraph and segments the time series into windows that can be used for consequence/impact analysis. For example, given a hypergraph time domain of $[t_0,t_f]$, a coincidence interval of $w$, and a shift $s$ (where $s\le w$), it is possible to create a set of windows ($W_k\in W$) with time intervals

As Donges et al. (2016) described, the ECA only considers the coincidence interval. The sliding window procedure adds more flexibility to defining coincidence windows using the shift $s$ parameter.

Each sliding window $W_k$ can be used to create a hypergraph snapshot $H_h$ by including any edge with a nonempty intersection between the sliding window interval $W_k$ and the edge’s interval set. The presence of two edges in a sliding window interval signifies independent concurrent occurrences. This is written aswhere $E_j$ = edge (i.e., hazard relationship) in the set of edges defined in the static hypergraph; and $I(E_j)$ = interval collection for edge $E_j$.

The collection of the hypergraph snapshots $\mathcal{H}=\{H_1,H_2,\dots ,H_t,H_n\}$, can be cast as a discrete dynamical process ($H_t⟼H_{t+1}$) to gain insight into the multihazard dynamics of the region in a considered timeframe (i.e., providing a temporal attribute). Each hypergraph snapshot can be mathematically expressed using a window incidence matrix (referred to in this paper as window incidence matrix to differentiate it from the static incidence matrix that is defined for a point in time) with elements $a_{i,j}$, where the column elements represent the edges $E_j$ and the rows represent the hazard sources $h_i$. The window incidence matrix depicts the occurrence of a primary or secondary hazard, linked to a particular single-hazard or multihazard relationship at a given window. Similar to a static incidence matrix, a window incidence matrix element is 0, if $h_i\notin Ej|W_k$, 1 if the hazard is primary, and $-1$ if secondary.

Stakeholders may choose to select a desired coincidence window (e.g., $w=1$ or 2 months) to characterize the joint occurrence of earthquake and flood hazards because it allows them to implement adequate flood-proof postearthquake recovery management strategies within the first two months of a big earthquake. The illustrative example presented in this paper demonstrates the significance of selected $w$ and $s$ on the temporal incidence matrix and capturing independent concurrent occurrences.

We consider a simple case where Fault 1 produces earthquakes with (moment) magnitudes $6\le \mathrm{M}\le 8$ at a constant rate of 0.01 per year. Fault 2 produces earthquakes with magnitudes $6.5\le \mathrm{M}\le 9$ at a rate of 0.025 per year. In other words, for the mainshocks, we assume for simplicity (and for illustrative purposes only) an HPP for interarrival time simulation and the doubly bounded exponential function (Cornell and Vanmarcke 1969) for earthquake magnitude simulation. Both faults are assumed to be strike-slip faults. Local seismic IMs—in terms of peak ground and spectral accelerations—are simulated according to the ground motion model of Boore et al. (2014) for illustrative purposes (any other relevant ground motion model for the region of interest could be used). The soil conditions are described later in this section. Note that, even if not very common, a strike-slip fault can generate a tsunami, as in the considered example (e.g., Ho et al. 2021). For instance, if a strike-slip fault is located near a shallow bay or a narrow inlet, the horizontal motion of the fault can push water into or out of the bay, creating a surge or a drawdown of the water level. This can amplify the wave height and cause flooding or erosion along the coast. An example is the 2018 Palu tsunami in Indonesia, triggered by a magnitude 7.5 strike-slip earthquake along the Palu-Koro fault (Ho et al. 2021). The earthquake caused a southward displacement of the seafloor along the Makassar Strait, which generated a long-wave tsunami that propagated into the Palu Bay and reached up to six meters high.

The aftershock sequence for the two faults is modeled using the BASS model (Turcotte et al. 2007). Other aftershock sequence modeling/simulation approaches [e.g., ETAS (Ogata 1998) and equivalent Poisson through Omori law (Shcherbakov et al. 2005)] can also be adopted here. Specifically, the BASS model is used to simulate the number, magnitude, occurrence time, and location of aftershocks, as shown in Eqs. (8)–(11), respectivelywhere $N_d$ = number of offspring earthquakes with magnitude $m$ greater than the specified minimum magnitude $m_{\min }$; $\mathrm{\Delta }m$ = difference between the magnitude of the parent earthquake $m_p$ and the largest aftershock [see modified Båth’s law (Shcherbakov et al. 2005)], $b$ = parameter of the Gutenberg-Richter (GR) distribution; $c$ and $p^{*}$ = temporal Omori parameters (Shcherbakov et al. 2004); $d$ and $q^{*}$ = spatial Omori parameters (Helmstetter and Sornette 2003); and $u_m$, $u_t$, and $u_d$ = cumulative distribution functions for the magnitude, time, and distance of the offspring earthquake, respectively. The GR-Båth, temporal Omori, and spatial Omori parameters used in carrying out BASS simulation in this study (see Table 1) are based on Turcotte et al. (2007). For the purpose of this study, similar parameters are assumed for both faults.

The considered liquefaction occurrence model combines the peak ground acceleration (PGA) at the site with the geotechnical parameters of the liquefiable site to simulate earthquake-induced liquefaction occurrence. The probabilistic distribution of PGA at the liquefiable site is estimated using the Boore et al. (2014) model, assuming $V_{s,30}$ equals $240\mathrm{m}/\mathrm{s}$. As shown in Fig. 2, the distance between the liquefiable site and Faults 1 and 2 are 10 km and 50 km, respectively. The assumed borehole data for the liquefiable site is shown in Fig. 6. The depth of the groundwater level at the site is taken as 1.5 m. Two soil types—clayey sand (SC) and poorly graded sand (SP) are assumed for the site. The soil density parameters are $14\mathrm{kN}/{\mathrm{m}}^3$ and $17.5\mathrm{kN}/{\mathrm{m}}^3$ for dry and wet soil, respectively. For this study, we assume the liquefaction occurrence time ≈ occurrence time of triggering earthquake.

The probability of liquefaction ($P_L$) can be estimated as (Idriss and Boulanger 2010)where $\mathrm{\Phi }$ = standard normal cumulative distribution function; and $F_s$ = factor of safety against the liquefaction for each layer. Other models (e.g., Cetin et al. 2004; Juang et al. 2012; Li et al. 2006) can be adopted as well.

The factor of safety is estimated aswhere ${\mathrm{CRR}}_{7.5}$ = cyclic resistance ratio calibrated for a M7.5 earthquake and is estimated using Eq. (14) (Idriss and Boulanger 2006); CSR = cyclic shear ratio and is estimated using Eq. (15) (Seed and Idriss 1971). The magnitude scaling factor (MSF) accounts for shaking duration impacts, and $K_{\sigma }$ is the overburden correction factor (Idriss and Boulanger 2006)where ${(N_1)}_{60cs}$ = equivalent clean-sand SPT blow estimated according to Idriss and Boulanger (2006)where ${\sigma }_{vc}$ = total vertical stress; ${\sigma }_{vc}^{\prime }$ = effective vertical stress at depth $z$; $a_{\max }$ = peak horizontal ground acceleration (i.e., PGA) at the ground surface; $g$ = gravitational acceleration; and $r_d$ = stress reduction factor that accounts for the flexibility of the soil column, estimated according to Idriss and Boulanger (2006).

A detailed probabilistic tsunami hazard analysis (PTHA) (e.g., Geist and Parsons 2006; Goda 2019; Horspool et al. 2014), which relies on numerical tsunami propagation methods, is required to simulate the tsunami runup height. Such analyses typically account for the topographic and hydrodynamic characteristics of the site and tsunami waves. For the intended illustration purpose of this section, rather than carry out a PTHA, we use the PTHA output from Goda (2019) [i.e., cumulative distribution functions (CDF) of tsunami wave height for different earthquake magnitudes] based on heterogeneous slip models for our tsunami severity model. Given that Goda (2019) provides CDFs for only eight magnitudes, we fit an exponential function to the median estimates of the CDFs for interpolation purposes [Fig. 7 and Eq. (16)]. As shown in Fig. 7, the tsunami wave height can be expressed as a function of earthquake magnitude only. We note that other studies (e.g., Abe 1995; Kulikov et al. 2005) have developed similar expressions using actual data from tsunamigenic earthquakes. A lognormal distribution with a standard deviation of 0.33 is associated with Eq. (16). Similar to the liquefaction hazard, we make a deterministic assumption that the arrival time of the tsunami wave ≈ occurrence time of the triggering earthquake (e.g., Carvajal et al. 2019)

In this example, climate change is assumed to influence the flood-inducing processes (i.e., snowmelt or extreme precipitation), resulting in a time-varying flood occurrence rate $\lambda (t)$ shown in Fig. 8. The hypothetical time-varying flood occurrence rate in Fig. 8 is informed by existing studies (e.g., Berghuijs et al. 2017; Mudelsee et al. 2004). In this paper, we make a simple assumption that the yearly maximum flood depth associated with climate-induced flood ($D_{\max ,1}$) follows a lognormal distribution with a time-varying relationship shown in Eq. (17) and a standard deviation of 0.4 m. Eq. (17) follows conclusions from studies (e.g., Kuok et al. 2022) that flood depth can increase by 1 m every 30 years. Similarly, other flood severity parameters (e.g., flood discharge) are not simulated for simplicity. Sophisticated hydrological and hydrodynamic models [see Schumann et al. (2013) and Yamazaki et al. (2013)] that account for topography and hydrodynamics can be used to simulate the severity and local intensity of flood hazards. However, such models may be computationally expensive. Simplified occurrence and severity models can be used, provided the trade-off between model simplicity, analysis time, and accuracy of results is justified.

Furthermore, we assume that any tsunami hazard with tsunami runup height $\ge 1\mathrm{m}$ triggers river overflow in the area. The maximum flood depth associated with the tsunami-induced flood ($D_{\max ,2}$) is hypothetically assumed to be dependent on the tsunami runup height using Eq. (18). We note that Eq. (18) is not based on actual data but only used here for illustrative purposes. Eq. (18) is assumed to provide a median estimate with a lognormal standard deviation of 0.5 m. It is also assumed that the arrival time of the tsunami-induced flood ≈ arrival time of the tsunami wave

The sliding window procedure can identify the joint occurrence of two independent events from different edges in a simulated event set. As mentioned earlier, the output of the sliding window procedure is temporal incidence matrices that highlight the joint occurrence of events over the specified sliding window. It is evident from Fig. 9 that the third primary flood event (month 533.12) occurs around a similar period as the seismic activity in Fault 1 (month 534.6). However, this could be missed in a numerical analysis without a sliding window analysis or if an unfavorably small window width is specified.

Let us assume decision-makers have defined two months as the upper bound limit to define the joint occurrence of two hazards for planning purposes. Fig. 10 presents the hypergraph snapshots generated from the sliding procedure with two combinations of $w$ and $s$ around the occurrence time of the third primary flood and seismic event in Fault 1. As shown in Fig. 10(a), by adopting $w=2$ months and $s=1\mathrm{month}$ [i.e., $W$ (in months) = $\{[0,2],[1,3],[2,4],\dots ,[598,600]\}$] we capture the joint occurrence of the third primary flood event and the seismic activity of Fault 1 in $H_{534}$ {i.e., window [533,535]} [see Eq. (19)]. Alternatively, if we adopt $w=1\mathrm{month}$ and $s=1\mathrm{month}$ [i.e., $W$ (in months) = $\{[0,1],[1,2],[2,3],\dots ,[599,600]\}$], the joint occurrence of both primary hazards is missed [see Eq. (20) and Fig. 10(b)]

Information derived from the sliding window analysis (i.e., hypergraph snapshots and window incidence matrix) can be used to design or assess multihazard management strategies. Furthermore, as mentioned earlier, the sliding analysis can help capture independent events occurring in quick succession. The window incidence matrices can be used, for instance, to define scenarios for post-disaster recovery analysis (e.g., Opabola and Galasso 2024) of civil infrastructure susceptible to multiple hazards.