Testing the Multivariate Hüsler–Reiss Model as a Practical Parametric Approach for Multiple River Flood Risk Assessment Using d4PDF Data: A Case Study in Kyushu Island, Japan

Due to a variety of hydrological extremes intensified in a changing climate, compound (multivariate) flooding events pose a serious threat across the globe. The escalating risk of compound flood is projected from a global scale (Zscheischler et al. 2018) to a national scale (Wahl et al. 2015) to a regional scale (Tanaka et al. 2020). Assessment of such multivariate extreme risks has been addressed using many multivariate models. In existing literature of water resources management, a number of dependence models have been applied and tested against multivariate data such as three variables (rainfall deficit, soil moisture deficit, and high temperature) for drought risk analysis (Sharma and Mujumdar 2022), high flow discharge of two tributary rivers (Neal et al. 2013), and two variables of compound flood risk factors (river flooding and storm surges) (Zellou and Rahali 2019). However, these studies targeted two or three variables due to the limitations of model flexibility and sample size (Kotz and Nadarajah 2000; Dey and Yan 2016). Some alternative approaches have been shown to overcome the former challenge. Heffernan and Tawn (2004) proposed a conditional approach to estimate the exceedance probability of multiple target variables as one exceeds a threshold. Owing to its flexibility, the nationwide assessment of concurrent flood risks in the United Kingdom (Keef et al. 2009), the European Union (Diederen et al. 2019), or the United States (Quinn et al. 2019) was implemented. Graphical models have also gained popularity for their flexibility in the dependence modeling of high-dimensional extremes (Engelke and Hitz 2020).

These approaches overcame the curse of dimensionality by avoiding from formulating multivariate extreme value distributions (MEVDs) in extremely high dimensions. On the other hand, from a family of MEVDs, a few practical models have been proposed and applied, trying to avoid the dimensional issues by constructing a model with fewer parameters. For example, Cooley et al. (2010) proposed a new construction of an angular density for the pairwise beta distribution, with only pairwise parameters and a global parameter ($2N+1$), and applied to spatial air quality data. They focused on the situation where threshold exceedance data were preferred over block maxima data. As introduced by Mostofi Zadeh et al. (2019), hydrological frequency analysis in practice prefers block maxima data because annual maximum data are a more commonly available form of data (Centre for Ecology and Hydrology 1999) and are away from the threshold selection issue (Solari and Losada 2012), which is also the case in the area of this study in Japan (Tanaka et al. 2002). However, such studies hardly show the explicit form of joint probabilities, which is more appropriate for multiple hazard risk assessment. Sharma and Mujumdar (2022) applied asymmetric logistic, Dirichlet, Hüsler–Reiss (HR), and pairwise beta models to meteorological and hydrological data (e.g., drought, monthly maximum temperature, rainfall deficit, and soil moisture deficit), showing that the HR model was flexible and outperformed the others in three-dimensional (3D) analysis. The HR model would be considered to take a central place of multivariate extreme value analysis, such as a multivariate normal distribution applied to general statistical problems, and moreover it is efficient, meaning that all the dependence parameters are obtained from the combination of pairwise parameters, while the target dimension is flexibly changeable, which are not satisfied in some popular models, such as the pairwise beta distribution (Husler and Reiss 1989). The graphical model proposed by Engelke and Hitz (2020) employs the HR model for building the dependence structure. The Heffernan and Tawn model is a conditional model, i.e., the aforementioned previous studies using this model calculated joint probability of all target locations in a certain region by multiplying conditional cumulative probability by the marginal cumulative probability of a single conditioned location as an origin. This framework avoids treating too many combinations. On the other hand, the multivariate HR model explicitly derives the whole function of multivariate cumulative probability, keeping simplicity by deriving dependence of all pairs from the bivariate HR model; therefore, this study explored the applicability of this model.

As a practical index of multiriver flood risk, instead of the joint probability of all rivers, this study examined the sum of subset joint probabilities (SSJP), defined as the probability that any combinations of $N$ rivers are flooded (exceeding such as 100-year design flow) over the target area (see the detail explanation and the formula of probability calculation in “Results”). For such a purpose, the parameter scalability in the HR model is useful because it can derive SSJP for various combinations only with pairwise parameters, while it is also a limitation in that, when any combination of pairwise parameters is inappropriate, the model cannot be constructed. Previous studies in related literature focused on the capability (flexibility or robustness) of spectral density (Boldi and Davison 2007; Engelke et al. 2015; Sharma and Mujumdar 2022). To demonstrate the aforementioned usefulness of the multivariate HR model more directly, this study estimated SSJP over a target region using the HR model.

In Japan, a heavy rainfall event during the Baiu (rainy) season in western Japan in 2018 and Typhoon Hagibis in 2019 caused multiple river flooding (MLIT 2018, b). However, such cases are rare in statistics, and thus hard to quantify as frequency based on observation data of river flow, which are available only for the past 30 years. The sample size issue has been addressed in the climate change research through an ensemble climate simulation approach. In parallel with the increased computing power and advanced climate sciences, several research institutes have successfully produced a large ensemble of climate simulation data. For example, in Japan, a database for policy decision-making for future climate change (d4PDF) has been jointly produced using the Earth Simulator, in coordination with various science programs [SOUSEI, TOUGOU, SI-CAT, and data integration analysis system (DIAS)] of the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) (Mizuta et al. 2017). This data set has been widely used in impact assessment of extreme weather events, such as river flooding, wind disaster, and coastal flooding (Ishii and Mori 2020). Kobayashi et al. (2020) developed a distributed rainfall-runoff model (1K-DHM) (Tanaka and Tachikawa 2015) for 109 major river basins (classified as Class A in terms of societal and economic importance by the Japanese government) and simulated annual maximum river discharges from d4PDF rainfall data with bias correction over Japan, which generated thousands of years of extreme river discharge data for periods of interest. Fully utilizing this large-ensemble, bias-corrected data set, this study aimed to overcome the sample size issue and clarify the asymptotic probability of simultaneously occurring flooding events of $l$ rivers among $n$ rivers. Tanaka et al. (2021) investigated a number of rivers with river discharges at the reference stations exceeding the design level, simulated river discharges from d4PDF rainfall data, and discussed geographical features via association rule analysis and graph theory. However, they neither analyzed asymptotic dependence nor discussed robustness of obtained dependence structure. Tanaka and Kitano (2022) analyzed the asymptotic dependence structure between two river basins using d4PDF, confirming that typical decadal record length is not enough to robustly identify tail dependency. On the other hand, its extension to multiple dimensions more than two including deriving SSJP has not been addressed yet.

The objective of this study was to propose a framework for estimating the multivariate probability of the number of flooded rivers (i.e., river discharge at the reference station exceeding the threshold return period) by deriving SSJP based on the HR model. Then, its applicability was tested in a case study in Japan using a large-ensemble river discharge data set of d4PDF. The country has been heavily damaged repeatedly. A heavy rainfall event during the Baiu season in western Japan in 2018 caused flooding in eight major rivers and resulted in 237 casualties, while Typhoon Hagibis in 2019 caused about 58 billion yen of insurance loss (MLIT 2018, 2019). In particular, we selected Kyushu Island where both severe rain bands in the Baiu season and typhoons often attack. “Data and Study Area” describes a large-ensemble climate simulation data set d4PDF and the study area. “Methods” introduces the methods starting with bivariate and multivariate HR model in “Bivariate Extreme Value Theory” and “Multivariate Hüsler–Reiss Model,” respectively, followed by the derivation of SSJP from various dimensional HR model “Developing the SSJP Model” and the simulation workflow in “Methodological Steps to Derive SSJP Based on Multivariate Hüsler–Reiss Model.” “Results” shows the fitting result of the SSJP with the HR model starting with a generalized extreme value distribution as a marginal distribution in “Performance of the Univariate Extreme Value Distribution,” bivariate dependence estimations and SSJP in “Two-River-Flood Probability” and the validation of derived SSJP in “Validation for Multiple River Flood Probability.” The characteristics of the obtained parameters and the applicability (limitations) of the framework are discussed in “Discussions.” Finally, “Conclusions” summarizes the discussions and limitations.

The data set of d4PDF was developed by a large ensemble of simulations of the past and future global climate (Mizuta et al. 2017) with the global climate model (AGCM-3.2S) at a spatial resolution of 60 km and downscaled with the 20-km nonhydrostatic regional climate model (NHRCM) (Sasaki et al. 2008) for Japan. The dynamically downscaled data in the past climate experiment contain 60-year climate simulations from 1951 to 2010 with 50 ensembles, amounting to 3,000 simulation years in total. As mentioned, observed river flow data in Japan are assumed to lack for multivariate extreme value analysis. Therefore, this study collected all ensembles of d4PDF for testing the multivariate HR model. Kobayashi et al. (2020) identified the annual maximum rainfall event in each year for 109 Class A river basins in terms of $N$-hour basin-average accumulation rainfall (BAAR), where $N$ is the design rainfall duration. They corrected model bias in BAAR from the observed BAAR using a quantile mapping method and simulated peak discharge from the bias-corrected rainfall data (defined as annual maximum floods) with a distributed rainfall-runoff model (1K-DHM) (Tanaka and Tachikawa 2015) for 109 major (named as Class A) river basins in Japan. As a result, the estimated design flow (corresponding to 100, 150, or 200 years in each basin) as the corresponding quantile from the simulated 3,000-year discharge was in a good agreement with actual design flow discharge set by the government [Fig. 10 in Kobayashi et al. (2020)]. This study employed the same 3,000-year data of simulated discharge in the past experiment in river basins in southern Japan and fit the MEVD framework to be described subsequently. This treatment was justified because the ensembles of boundary conditions of sea surface temperature and sea ice were all based on the observed data from COBE-SST2, and added perturbation was small compared with the observed value itself.

Simultaneous flood risks were investigated for 20 Class A river basins over the Kyushu region in Japan (Fig. 1). The hydrological summary of the target river basins is given in Table 1. The Kyushu region is an island isolated from the mainland of Japan. The island is composed of relatively small catchments (from 87 to $\mathrm{2,860}{\mathrm{km}}^2$) and often experiences heavy rainfall. In particular, the west side is heavily affected by seasonal frontal rain during May–July (known as Baiu), while the east side receives huge rainfall by typhoons during August–October. Since 2010, heavy rainfall events have caused serious flood disasters frequently, especially in the west side by frontal rain events (2012, 2017, 2018). Therefore, high dependence among these river basins is expected in Kyushu Island.

The performance of the constructed model was tested against the d4PDF data (annual maximum river discharges). Making the best use of d4PDF with large ensembles (3,000 annual maximum data), various block sizes were tested. First, the convergence of univariate extreme values to GEV was explored. The river discharges in all combinations of two river basins were modeled by the bivariate HR model, which was validated in terms of the CDF of the angular component. Using the estimated dependence parameters, multivariate HR models and the resulting SSJP among multiple river floods were estimated and compared with the empirical CDF of SSJP obtained from the d4PDF data. In all analyses, R package evd (Stephenson 2002) was used. Parameters were fitted with the maximum likelihood method.

In the preceding validation, an empirical CDF was created by counting the number of river basins where annual maximum river discharge exceeded the threshold return periods (100, 200, 500, and 1,000) in d4PDF data. Although the developed SSJP framework based on multivariate HR model was not fit to this higher-dimensional ($\ge 3$) empirical CDF and estimated by combining only pairwise dependencies in the bivariate HR models, the result showed acceptable agreement. This suggests that the proposed framework fairly works for estimating flood risk over catchment.

Whether or not the multivariate HR model is constructable from bivariate HR models depends on the combination of bivariate dependence parameters. Given the constraint of partial correlation ${\rho }_{k;i,j}$, the dependence parameters needed to satisfy the following conditions:

The preceding parameter space is shown as a contour of ${\rho }_{k;i,j}$ on the 2D coordinate between $\log (q_{ik}/q_{ij})$ and $\log (q_{jk}/q_{ij})$. According to Fig. 10, the conditions for Eq. (17) indicated that if one parameter was relatively small or large, the others must be similar. The similar value of the three inferred previously was located at around the origin in the figure. The higher the dimension, the harder it was to satisfy Eq. (17) for all the combinations of the three dependence parameters. Thus, the applicable number of rivers for the HR model declined, though the number actually increased from three to five dimensions simply due to the increasing number of total combinations. Eq. (17) is a simple but useful index to determine whether data satisfied the minimum requirements of the HR model.

The histograms of mean and coefficient of variation (COV) of $\chi$ for the modeled combinations and others (some partial correlations were negative) in the range of three to seven dimensions are shown in Figs. 11 and 12, respectively. Overall, the mean values of $\chi$ were small. The modeled combinations showed relatively higher correlation coefficients. The COV values of $\chi$ were smaller for the modeled combinations. This tendency was clearer in larger dimensions. As discussed previously, a combination of similar partial correlations (around the origin of Fig. 10) falls within the parameter space for the HR model (the outer line in Fig. 10). On the other hand, among the applicable combinations in Figs. 7 and 8, all the river basins were not necessarily from the same major group in Fig. 5, but a few were from the other (e.g., one example per river basin was River 91, 93, 97, 98, 99, 101, and 105). In such a case, a partial correlation in some combinations fell into the limits of the allowable ranges in Fig. 10. In case of nested logistic (Coles and Tawn 1991) or nested negative logistic models (Kitano 2022), they can only express the dependence conditions where many pairs among all bivariate combinations share the same dependence structure, which is somehow a relaxed model rather than a logistic and negative logistic model with a single parameter. The HR model is considered as a more flexible model, whose total number of parameters is the number of pairs, and the relation of parameters values are less restrictive while still simple enough to analytically derive SSJP. In other words, the HR model can be an extension model by relaxing the nested model.

According to Fig. 5, the target catchments in Kyushu Island was classified into two groups of mutually correlated river basins: one faces the East China Sea (west side, Rivers 90–101) and the other is located on the Philippines Sea side (east side, Rivers 102–107). Heavy rainfall on the east side of the island is typically supplied by typhoons that generated on a path through the Philippines Sea, while that on the west side is more characterized by a seasonal frontal rain band (called Baiu) from June to July (MLIT 2015). The proportion of occurrence season of extreme river discharge ($>100\mathrm{years}$) in d4PDF between the Baiu season (June or July) and the typhoon season (August or September) was investigated in Table 3. Among 20 river basins, Rivers 90–101 had more events in the Baiu season, while the typhoon season is more dominant for Rivers 102–107. Rivers 108 and 109 show a relatively higher ratio of occurrence in the Baiu season; therefore, the correlation $\chi$ was lower with Rivers 102–105 in the same group. The Ministry of Land, Infrastructure, Transport and Tourism (MLIT) reported that River 108 experiences heavy rainfall in both seasons (MLIT 2007b), while flood events in River 109 have been historically caused by typhoons (MLIT 2007a). Although inconsistency of flood occurrence season in River 109 may be caused by systematic bias of typhoons in d4PDF, overall characteristics of extreme flood seasons in Kyushu Island are expressed in d4PDF, which are clearly shown as mutually correlated groups of river basins in Fig. 5. This classification agreed with the rainfall-based hydrological regionalization by MLIT (2015). Hydrological classification has been also utilized in other case studies to efficiently construct high-dimensional dependence models in a literature (e.g., Quinn et al. 2019). As shown in Fig. 10, it is hard for the combination of high- and low-dependence parameters to satisfy the constraint of partial correlation coefficients. Therefore, it is suggested that a priori hydrological classification of target rivers may improve the feasibility of deriving SSJP using the multivariate HR model.

Simulated river discharge data from d4PDF are known not to converge to the extreme value theory with the block size shorter than 20 years (Kitano 2022). Therefore, the uncertainty of the derived SSJP may be wider when applied to past observation data. Another limitation is that as the number of target river basins increases, the multivariate HR model gets harder to construct. As a result, a nonnegligible number of trials were discarded, when CIs could not be created from the rest of the samples because it leads to biased estimation. It is clarified that the strong regulation of partial correlation makes it hard to derive CIs when the HR model is applied to derive SSJP. Still, CIs estimated from empirical SSJP added on Figs. 7 and 8 provide supportive information of sampling uncertainty.

Intensified by global climate change, the compound risk of the co-occurrence of hydrological extremes, such as multiple-river flooding or simultaneous river and coastal flooding, receives high public attention. This study devised a multivariate HR model to estimate compound exceedance probabilities as the SSJP for a large number of river basins and tested its applicability based on extreme river discharge data from a large ensemble ($>\mathrm{1,000}\mathrm{years}$) of climate simulation data (d4PDF). The obtained findings are summarized as follows:

The successfully constructed multivariate HR model was composed of the river basin pairs with similar dependence structure. Nevertheless, the HR model was found to allow a different dependence structure to some extent and thus was more flexible than the nested logistic or negative logistic model, which requires the same dependence structure among variables.

Overall, this study demonstrated promising results for applying the HR model to the co-occurrence of multiple river flooding under certain conditions. The large ensemble of d4PDF substantially ensured the convergence of the data to asymptotic theory, which contributed to robust compound risk assessments. In a bivariate case, Kitano (2022) has already revealed that observation record length is not enough to robustly identify tail dependency. Owing to the sample size secured by d4PDF, the preceding applicability (and limitation) of the SSJP based on the HR model was derived. The systematic bias in d4PDF was removed by Kobayashi et al. (2020); however, estimated tail dependency is based on climate simulations. Future study will cross-check the joint flood risk based on both observed and simulated data and explore climate change impact by applying the methodology to future climate projections.

Some or all data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request.

This work was supported by the MEXT Program for the Advanced Studies of Climate Change Projection (SENTAN) Grant No. JPMXD0722678534 funded by the Ministry of Education, Culture, Sports, Science, and Technology, Japan (MEXT). This study used d4PDF produced jointly with the Earth Simulator by science programs (SOUSEI, TOUGOU, SI-CAT, DIAS) of MEXT, Japan. This data set was collected and provided by the Data Integration and Analysis System (DIAS) (Project No. JPMXD0716808999), developed and operated by MEXT.