Non-Monotonic Relationships between Return Periods of Precipitation Surface Hazard Intensity

: Hazardous surface processes such as ﬂoods and mass movements are often induced by a common trigger such as extreme precipitation. The relationship between the intensity of the trigger and the surface hazard is generally assumed to be monotonically increasing (increasing precipitation never decreases hazard intensity). The validity of this assumption of complex multi-hazard events has not been thoroughly investigated. In this research, the relationship between cumulative precipitation and hazard intensity was investigated by a simulation of 50 return period precipitation events on the Carribean island Dominica. Here, several tropical hurricanes have induced events with (ﬂash) ﬂoods, slope failure, debris ﬂows and landslides within the past decades. Results show that complex multi-hazard interactions break the common assumption for the relationship between trigger and hazard intensity. In particular, landslide dam formation and mass movement dilution result in hazard intensities that are not a one-to-one increasing function of trigger intensity. Spatial variability in this behavior is quantiﬁed using a rank-order correlation coefﬁcient between trigger return period and hazard return period. Since trigger and hazard return periods are, in the study case, not approximately equal, the hazard for a speciﬁc location can not be classiﬁed based on trigger return period. This has implications for risk calculation and decision making related to disaster risk reduction.


Introduction
Storm-related precipitation can induce hazardous surface processes such as flash flooding and riverine flooding. The intensity and extent of the surface hazard (e.g., flood) dynamics depend directly on the spatio-temporal distribution of the precipitation [1,2]. For this reason, the assessment of the hazard from flood processes requires quantification of storm properties [3]. The assessment of hazard and risk for storm-induced surface hazards generally consists of three major steps [4]. First, analysis of past storm events must provide a relationship between the return period and the storm properties such as cumulative precipitation. Secondly, surface process dynamics are estimated by physically based modeling of flood flow. Finally, using the obtained return period (or yearly probability of occurrence) and hazard intensity, additional calculations for impact and risk are carried out.
A major assumption in these processes is the monotonically increasing (never decreasing, monotonic) relationship between the intensity of the storm and the intensity of the surface hazard. It is assumed that increasing precipitation results in increased flood intensity (e.g., flow height). This relationship is both theoretically and experimentally confirmed for numerous surface hazards [5][6][7]. In the case of floods and storm precipitation, direct increased in rainfall results in higher runoff and increased flooding, as is confirmed by theory and observations. For co-seismic landslide initiation, a similar relationship holds. Generally, an increasing peak ground acceleration results in higher co-seismic landslide density [8]. Finally, in the case of hurricanes, categorization generally takes place on the basis of maximum wind speeds [9], which is assumed to directly relate to increased impact [10]. This monotonic relationship between storm intensity and the intensity of hazardous surface processes such as flooding is used in a variety of calculations related to hazard and risk assessment. First, risk is generally conceptualized as the combination of hazard, exposure and vulnerability [11]. The calculations of total risk therefore involves the definition of the return period P(T), which is then used to calculate the integral over all possible events (R = ∞ 0 P(T)V(I(T))dT, where V is the vulnerability and I is the intensity). Secondly, the categorization of natural disasters is often performed by means of finding the return period for the event. The assigned category can then influence policy, measures and insurance payout. Finally, prediction of future hazard and impact requires the design of scenarios for these events [12]. For storms, this commonly occurs by means of linking return periods of interest to design storms. The predicted hazard then influences decision making related to disaster risk reduction and spatial planning.
Despite the importance of the monotonic relationship between storm intensity and surface hazard intensity, the extension towards multi-hazard events has not been properly investigated. Extreme storms can, besides flash floods and riverine flooding, trigger slope instability and mass movement runout. In addition, these processes might be triggered simultaneously and interact, resulting in new types of dynamics and hazard types, such as debris floods and landslide dam formation [13]. Such multi-hazard events are generally recognized to increase total hazard beyond their individual contributions [14]. An example of these types of multi-hazard events can be found in the impact of Hurricane Maria on Dominica [15]. For examples of such multi-hazard models, see r.avaflow [16] or OpenLISEM Hazard [13]. In numerous studies, the complex non-linear influence of multi-hazard behavior on spatial hazard intensities are noted [17]. While research on the dynamics of multi-hazard events has resulted in integrated multi-hazard modeling tools, the influence of multi-hazard interactions on the relationship between trigger intensity and surface hazard intensity has not been properly studied.
In this research, we aimed to quantify the potential non-monotonic complexities in storm-induced multi-hazard events featuring floods, debris flows and landslides. This is conducted by means of the physically based multi-hazard simulation of a large set of return period events. The studied event is the impact from hurricane Maria on Dominica, which featured various types of mass movements and (flash) floods. From the results, the relationship between trigger and hazard intensity was studied. Additionally, we quantify the spatial variability in this relationship and the primary causes within our studied set of events. In Section 2, we describe the study site in more detail. In Section 3, the multi-hazard modeling tool is explained. Finally, Sections 4 and 5 provide modeling results, analysis and discussion.

Study Site
The study site is the South-Eastern Grand-Bay catchment within the Caribbean island of Dominica ( Figure 1). Dominica, an island of volcanic origin, features steep terrain across its 750 km 2 surface area. The volcanic peaks reach over 1400 m in altitude. The island consists mainly of volcanic deposits but features some small zones of uplifted limestone. On top of the volcanic tuff and ignimbrite rocks clay-like weathering products form a shallow soil layer. Throughout the island, variation in soil properties depend on the age and weathering stage of the soils. The surface of the island is covered predominantly with dense tropical forests and minor areas of agroforestry and urban development.
Within the past decades, several major tropical hurricanes have impacted Dominica. Hurricane David (1972), Erika (2014) and Maria (2017) have impacted the Island directly leading to a variety of surface processes including (flash) flooding, slope failures, landslides, debris flows and mass erosion. Erika and Maria lead to a cumulative precipitation of over 800 and 500 mm, respectively, within 24 h. For Maria, hazardous flow processes were present throughout the island, with over 20,000 mass movements mapped, and wind impacted over 90 percent of the roofs on the island. A variety of multi-hazard interactions occurred throughout the Island for all of the mentioned hurricanes [18]. First, the occurrence of mass release from slope and flash flooding results in mass movement dilution, which can prolong runout significantly. Secondly, the combining of solid and fluid flow volumes results altered flow dynamics. In particular, there is increased drag, pressure, altered erosion and the occurrence of other interfacial forces within the resulting two-phase flow. Finally, reports and simulations indicate partial or full landslide dam formation occurred in numerous river-branches.

Design Events
For the study site, a series of design storms was created with return periods between 1 and 50 years. For each of these, the multi-hazard impact was simulated with an identical landscape state. The derivation of the design storms was part of the CHARIM Caribbean multi-hazard analysis project [19]. This started with analysis of a 20-year rainfall record from St. Lucia. Daily cumulative precipitation was assumed to represent the cumulative precipitation of distinct events. A rank-ordered General Extreme Value analysis (GEV) on yearly maximum event rainfall resulted in a relationship between return period and maximum precipitation. The GEV curve and 20-year return period event rainfall is shown in Figure 2. For each return period, a temporal precipitation graph is required. This timeseries should represent the characteristics common for the area. Since the aim of this research is to isolate the influence of multi-hazard interactions, the rainfall curves are chosen to have minimal variability. For a 20 year return period event, a minute interval rainfall intensity curve was available from a distribution fitted to an actual event. This curve was then re-scaled to represent the precipitation for the other return periods. By removing the temporal changes in the design storms, a simplification is made. The design event does not necessarily represent a best estimate of the storms; however, this simplification allows for better isolation of the influences of multi-hazard interactions on the relationship between hazard intensity and storm precipitation. Variations in the temporal and spatial patterns of precipitation have been known to influence flood intensity significantly and this issue has been studied intensively by others [20].

Modeling Setup
Simulation of the surface processes are carried out using the OpenLISEM MultiHazard model. This open-source tool provides simulation of catchment-scale hydrology, flow, (flash) floods, slope stability and mass movement hazards. A short overview of the model is provided in Figure 3. Hydrology was implemented according to the OpenLISEM catchment hydrology description [21]. This includes a two-layer Green and Ampt infiltration model combined with Darcy depth-averaged ground water flow. A generalized two-phase fluid/solid flow description was used, which was adapted from Pusadaini [22] to implement shallow surface flow [13]. Slope failures were simulated based on a real-time interaction with hydrology using the iterative failure method [13]. Finally, erosive processes for mass movements were implemented using Takahashi's set of entrainment equations [23].
Model input data consist of 21 inputs in total [13]. All of these are derived automatically from several basic input maps. An overview of the spatial and non-spatial model inputs is provided in Table 1  Analysis of the hazard intensities was performed by means of calculating the total momentum p = h f v f ρ f + h s v s ρ s with p the momentum, h the height, v the velocity and ρ the density of the material. The subscripts f and s indicate the fluid and solids phases of the flow, respectively. Total momentum provides a robust indicator of total hazard intensity for its inclusion of multiple physical properties of the flow. This is contrast to single-parameter intensity measures such as flow height or flow velocity. All of these types of intensity indicators are used in the vulnerability and impact assessment; however, calculations performed for multi-phase flow processes typically focus on multi-parameter indicators such as momentum or impact pressure [25]. To better isolate the processes influencing the relationship between trigger and hazard intensity, the analysis of flow height as intensity is similarly shown in the result section.

Results
Within this section, the general multi-hazard simulation results are shortly described, by using the 20-year return period simulation as an example. Afterwards, the relationship between trigger intensity and hazard intensity is analyzed by means of (i) surface hazard return period maps; (ii) spatial quantification of rank-order correlation between trigger and hazard return period; (iii) the distribution of the rank-order correlation values; (iv) plots of the relationship between hazard and trigger intensity for several locations.

Simulation of the Multi-Hazard Event
As a reference indication of the general interacting behavior of multiple hazardous land surface processes as simulated based on high-intensity rainfall input, the 20-year return period simulation output is shown in Figure 4. For calibration and validation, an analysis of comparison to mapped surface processes is provided by (Bout et al., 2020).
The basic simulation output is a rough reconstruction of the impact of Hurricane Maria, as the landscape description used as model input is derived from pre-Maria (2017) data. Flash flooding occurs in every major branch of the hydraulic network in the area, resulting in major flooding near the coastal outlet. Within the simulation, the slopes throughout the area are destabilized by means of infiltration and pore pressure increase, resulting in numerous mass movements. Field observations of shallow soil-slips in highly over-saturated soils support this cause for slope instability. The mass movements are distinct in type based on the material properties, terrain and failure depth. A part of the mass movements behave landslide-like in terms of retaining of some structure and limited spreading in its movement. The vast majority behaves debris-flow-like due to the high fluid content and steep slopes. The debris flows predominantly lose momentum once reaching the channel network, where slopes are reduced compared to the generally steeply sloped terrain. The strong fluid flows can, for the majority of locations, provide sufficient shear stress to drag or entrain material. Finally, the majority of deposits end up within the main river branch or near the outlet. Three major types of interactions occur within the simulation. The first, as mentioned above, is the dilution of mass movements by surface water flow. The additional water influences basal inter-granular pressures and increases lateral spread in the flow equations. The second interaction is the drag and entrainment of material. With sufficient water flow, soil or deposited solids can be moved together with the fluid. This can increase flow momentum as average density increases and friction forces behave differently within the model. The third interaction is the blocking of river branches by means of landslide dam formation. This occurs in various locations, most notably near the town of Pichelin, in the river segment that encloses the town on the West side.

Relationship between Trigger and Hazard Intensity
With the variety of processing taking place within the simulations, a total of 50 simulations with return period between 1 and 50 years have been carried out. For some representative types of locations, the hazard intensities as a function of the precipitation return period are graphed in Figures 4 and 5. A first and immediate observation within Figure 5 is the non-linear complexity in the relationship between trigger (precipitation) intensity and hazard intensity (flow momentum). While traditional flood and flow processes are expected to increase, potentially non-linearly, in a one-to-one relationship; this is not observed uniformly with the introduction of multi-hazard interactions. Locations three and four follow this assumption. These locations are located within the main river and are relatively uninfluenced by local multi-hazard interactions. Instead, intensity here is the product of the combined output of all upstream subcatchments. Because of this, the complex variability averages and only the trend of increasing flow volume is visible. A second, more common phenomena, are threshold effects leading to a large increase in hazard when a specific threshold trigger intensity has been passed. This can occur, for example, in flood modeling when an area is protected from flooding by an obstacle of a specific height. Once the water reaches this height during a simulation, the protected area can experience a rapid increase in flood height. This effect is observed within the simulation for several small areas. An additional strong threshold effect is the initiation of slope failures. These do not occur if precipitation is insufficient to reduce the factor of safety below one (see locations 9 and 10). The majority of slope failures therefore only occur beyond a return period between 2 and 7 years.
Five types of dynamics are observed within the relationships between triggered intensity and hazard intensity: (i) predominantly increasing (locations 3 and 4), (ii) approximately constant (locations 7 and 8), (iii) predominantly decreasing (locations 5 and 6), (iv) threshold increasing (locations 9 and 10), (v) no clear relationship (locations 1 and 2). Examples of these curves are shown in Figure 5. The predominant causes for these types of dynamics are listed in Table 2. Table 2. Observed types of relationship between precipitation intensity and surface hazard intensity.

Relationship
Observed Causes

predominantly increasing
The expected behavior for areas predominantly impacted by (flash) floods. As increase in water results in increased runoff, the converging water that forms the flash floods is similarly increased.
approximately constant Observed in particular for areas impacted with landslide type-movements. With increased trigger intensity, there is no strong influence on landslide flow momentum. Generally, slope failures reach the interface between soils and bed-rock, and additional precipitation does not increase the released volume. Instead, increased trigger intensity influences the spatial extent of impact by increasing the number of slope failures.
predominantly decreasing This type of relationship between hazard intensity and trigger intensity is found in case of shallow slope failures that move down on a diverging (laterally convex) slope. Here, the increase in trigger intensity (increased precipitation) results in higher water flow. This dilutes the mass movement resulting in spreading of the flow as fluid-pressure forces become more important. The spreading reduces the flow height and momentum for a particular location, although it increases the total exposed area.

threshold increasing
This type of relationship is found in cases where a threshold effect is significant, such as found in case of slope failures or areas protected from floods by a barrier. Additionally, the breaching of landslide dams shows threshold behavior as material is entrained only above a specific shear stress, which depends on flow heights, velocities and densities.
no clear relationship This is predominantly found for locations with interactions that depend on timing of individual processes such as landslide damming of rivers, breaching and entrainment. The potential damming of rivers depends on the arrival time of flash flood waves and deposition behavior of the mass movement. With sufficient water flow, breaching occurs, which prevents formation of a small reservoir behind the blocking. These time-dependent interactions do not show clear patterns but alter the relationship between trigger intensity and hazard intensity in complex ways. Figure 6 shows a spatial indication of the relationship between the return period of the triggering precipitation and the surface hazard intensity.

Spatial Variability
There are locations for which the return period for the hazard intensity is identical to the trigger return period. Despite this, the difference between these two values predominantly lies within several years. In particular, slope failure and mass movement impacted areas show large deviations. Downstream areas deviate less as influences from upstream multi-hazard interactions are relatively less important when outflow from multiple area join.

Quantification of Relation
A quantification of the deviation from a one-to-one monotonically increasing relationship between trigger intensity and hazard intensity can be made using a rank-order correlation coefficient. Spearman's rank-order correlation coefficient is here used to show the spatial variability Figure 7. The spatial patterns in the value of Spearman's rank-order correlation coefficient for the return period of the trigger and hazard intensity indicate where the multi-hazard interactions break the traditional assumption of one-to-one increasing hazard. A value of 1 indicates this assumption is perfectly held for all simulation. A value of −1 instead indicates a monotonically decreasing relationship. A value of 0 indicates no clear relationship is found. For the rank-order correlation, the absolute values do not alter the value of the coefficient, only the classified ranks are used in the calculation. Thus, as long as each consecutive hazard intensity increases, the shape of the relationship can be ignored and the coefficient is 1. The results show that slope failures, mass movements and complex terrain are the predominant causes of complexity. On the other hand, influences in downstream are relatively minor. This can be well explained by the effective averaging of the influence from each upstream sub-catchment.

Influence on Impact
The average rank-order correlation between trigger intensity and hazard intensity for several land use/cover or terrain types is shown in Table 3. Complexities in the relationship between precipitation and surface hazard intensity occur more strongly in the categories steep or urban surfaces. The predominant reason for this can be found in the influence of mass movements on the rank-order correlation between precipitation and surface hazard return period. The influence on built-up areas are explained by the building locations. Within the Grande Bay area, buildings are often located near the main channels, as the steep terrain prevents construction in other areas. These areas, such as the town of Pichelin, are influenced by multi-hazard interactions such as mass movements and landslide dam formation in relatively higher frequency compared to other areas. Finally, buildings are generally not located in flat areas prone to frequent flooding. The frequent hurricane-induced surface inundation prevents reasonable livelihood in these locations.

Influence on Risk Calculations
The assessment of hazard and risk is an important tool in decision making for disaster risk reduction. In the calculation of risk, linking hazard intensities to return periods is a crucial step. The formulation of risk as hazard times exposure times vulnerability is still locally valid (R = ∞ 0 P(T)V(I(T))dT). Additionally, regional calculations of total risk remain valid depending on the manner in which the integral over all possible events is approximated. When multiplying event probability with event damage, accuracy depends on the discretization of the return periods; however, the complex relationship between trigger and hazard intensity might stress the importance of using a large number of return period scenarios in risk calculation. For event-specific risk calculation, the return period of the trigger no longer indicates the return period of the risk value. This results from possible deviations of the return period of local impacts and damages from the return period of the trigger.

Influence on Event Categorization
In the presented study case, the classification of events based on the trigger return period is not valid. An event with low return period trigger might have a high return period value in specific locations and low in others. This spatial variability complicates a simple classification. In practice, the return period classification of events are often used in decision making and insurance [26]. In the case of hurricanes, return periods can determine both both policy in post-disaster reconstruction and insurance payout [27,28]. If the event is multi-hazard, an additional complexity arises from the classification in types. In the presented study case, debris flows, landslides and floods are not easily separable. Within the literature, there has been similar critique on the simplified usage of hurricane categories that are used for insurance payout calculations [29]. In particular, hurricane categories commonly depend on a single parameter and do not reflect real impact [10,30]. While these simplifications can be useful when simple rule sets are required for fast and transparent decision making, the complexity of multi-hazard intensity and trigger intensity shows that such simplifications can result in significant errors. If and how such complexity might be included in fields such as insurance, policy and urban planning must be further investigated.

Limitations of the Used Method
The multi-hazard modeling setup used in this work has been applied and validated on several study sites [13,31]. Still, several limitations exist because of the multiple hazards included in the model. Uncertainty is increased by the increase in input parameters and processes. The progression of uncertainty through these types of models remains to be investigated further; however, the results of this work and others show a significant spread in outcomes not typically observed in models for individual processes [17]. Usage of model results for hazard assessment and application might therefore require analysis of probabilistic modeling ensembles [32]. Another important limitation is the lack of additional processes such as wind, windthrow of trees, the influence of tree debris on flow dynamics and storm surging. All of these processes either influenced the impact on Dominica during tropical cyclone Maria, or interacted with other hazards to alter their dynamics [33]. Finally, the return period analysis of precipitation events was carried out using General-Extreme-Value theory, which has been successfully applied in many cases, but might be limited when applied to Hurricane systems [34].

Conclusions
The relationship between trigger and surface hazard intensity is commonly assumed to be one-to-one monotonically increasing. In this work, an ensemble of simulations with rainfall of increasing return period for the tropical Island Dominica has shown that multihazard interactions break this assumption severely in local areas. In particular, slope failure and mass movement impacted areas, including river segments prone to landslide dam and reservoir formation deviate in their relationship between trigger return period and hazard intensity return period.
A major influence of the found complexities is in the usage of trigger return periods for hazard classification. These classification are often used in decision making and insurance applications. There, considerations of local variability in hazard intensity return period can significantly alter the validity in existing methodologies; however, further research is required to investigate the extent to which the findings of this work impact other types of hazardous processes and areas.