# A Bayesian Framework for the Multiscale Assessment of Storm Severity and Related Uncertainties

Juliette Blanchet, Victor Mélèse

► To cite this version:

Juliette Blanchet, Victor Mélèse. A Bayesian Framework for the Multiscale Assessment of Storm Severity and Related Uncertainties. Journal of Hydrometeorology, 2020, 21 (1), pp.109-122. 10.1175/JHM-D-18-0254.1. hal-03087660

HAL Id: hal-03087660

https://hal.science/hal-03087660

Submitted on 12 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Bayesian Framework for the Multiscale Assessment of Storm Severity and Related Uncertainties

JULIETTE BLANCHET AND VICTOR MÉLÈSE

Université Grenoble Alpes, CNRS, IGE, Grenoble, France

(Manuscript received 20 December 2018, in final form 18 November 2019)

## ABSTRACT

This article proposes a statistical framework for assessing the multiscale severity of a given storm at a given location. By severity we refer to the rareness of the storm event, as measured by the return period. Rather than focusing on predetermined spatiotemporal scales, we consider a model assessing the return period of a storm event observed across the continuum of durations and areas around a focus location. We develop a Bayesian intensity-duration-area-frequency model based on extreme value distribution and space-time scale invariance hypotheses. The model allows us to derive an analytical expression of the return period for any duration and area, while the Bayesian framework allows us by construction to assess the related uncertainties. We apply this framework to high-resolution radar-rain gauge reanalysis data covering a mountainous region of southern France during the autumns 2008–15 and comprising 50 rain events. We estimate the model at two grid points located a few kilometers apart on either side of the mountain crest, considering spatiotemporal scales ranging over 3–48 h and 1–2025 km². We show that at all scales and for all significant events, the return period uncertainties are skewed to the right, evidencing the need of considering uncertainty to avoid systematic risk underestimation. We also reveal the large variability of the storm severity both at short distance and across scales, due to both the natural variability of rainfall and the mask effect induced by the mountain crest.

## 1. Introduction

Mediterranean storms triggering floods have long been a source of heavy damage caused to people and infrastructure. What is the frequency of occurrence—or equivalently the return period—of such devastating storms? For rain falling nonuniformly in space and time, the return period of a given storm depends on the location, the duration, and the area over which it is observed. This article deals with these three sources of variability.

Extremes are rare by definition, so the computation of extraordinary return periods—which are of most interest—faces a problem of lack of data. Different statistical methods have been proposed to strengthen the robustness of return period estimates based on short time series. Among them, the regional frequency analysis (RFA) is based on the identification of homogeneous subregions in terms of rainfall distribution, allowing the pool together of extreme rainfall data of a given subregion, and thus to base return level estimation on longer series.

Norbiato et al. (2007) used RFA to compute the return period of point rainfall during a storm that occurred in August 2003 in the eastern Italian Alps. The authors illustrated in particular that the rareness of point rainfall (i.e., the return period) depends on the considered accumulation duration, due to rainfall temporal intermittency. Liu et al. (2015) also used RFA to study a storm that occurred in February 2007 in Jakarta, Indonesia. They highlighted in particular the strong spatial variability of point rainfall rareness induced by the spatial intermittency.

Another way to strengthen the robustness of return period estimates is to account for the structural properties of rain fields, in order to base estimation on a pool of extreme rainfall depths observed at different scales. Several articles applied this idea for point rainfall considered at different durations. For example, Bougadis and Adamowski (2006), Van de Vyver and Demarée (2010), and Mélèse et al. (2018) used the temporal scaling properties of rainfall to derive analytical expressions of intensity-duration-frequency (IDF) curves representing the relationships between rainfall intensity, duration and return period of an event. The studies accounting for both the temporal and spatial scale—that is,

<footnote>*Corresponding author:* Juliette Blanchet, juliette.blanchet@univ-grenoble-alpes.fr</footnote>
dealing with extreme values of areal rainfall consid-
ered over different durations and areas—are much more
scarce, partly due to the lack of spatial rainfall measure-
ments. One of these is the study of De Michele et al.
(2001), which used the space–time structural prop-
erties of the rain fields to build analytical expres-
sions of intensity–duration–area–frequency (IDAF)
relationships in the Reno basin, Italy. Mélèse et al.
(2019) proposed for the region studied in this article a
new IDAF model that can be seen as a compromise
between the structure-based model of De Michele
et al. (2001) and a purely data-driven model, such as
the framework proposed in Overeem et al. (2010).
Additionally to providing robust return levels esti-
mates, the IDAF relationships allow to compute the
return period of an event for any duration and any
area in the applicable range of scales. Thus they allow
to assess the multiscale rareness of a given storm at a
given location. Following Ramos et al. (2005), we refer
to the rareness of a storm observed over a continuum
of spatiotemporal scales as the storm “severity” (or
multiscale severity). Ramos et al. (2005) considered the
storm severity in a regional way, by focusing on the
return level of the maximum rainfall value over a given
region. They used the IDAF relationships to study the
severity of three storms that occurred in a region of
300-km² area, around the city of Marseilles, in south-
eastern France. They revealed in particular that areal
return period estimates can vary from ordinary values
(1–2 years) to exceptional values (100 years) depending
on the considered durations and areas. Also, for a given
duration, both ordinary and exceptional point return
period estimates can coexist a few kilometers apart due to
spatial intermittency. Ceresetti et al. (2012) extended the
work of Ramos et al. (2005) to a 3000-km² region located
in southern France, using the IDAF model developed in
De Michele et al. (2001). By considering the storm se-
verity in a regional way, neither Ramos et al. (2005) nor
Ceresetti et al. (2012) could assess the spatial variability
of the severity across scales.

Another issue when dealing with return period esti-
mation is that uncertainty is unavoidable both because
models are imperfect and because data have a limited
length. Uncertainty is particularly large for the largest
return periods, which in turn are usually of most interest.
Evaluating uncertainty is required to avoid a system-
atic underestimation of the probability of occurrence of
disastrous events (Coles and Pericchi 2003). To the best
of our knowledge, the only study providing uncer-
tainty estimates of IDAF relationships is Overeem et al.
(2010), where uncertainty is estimated using bootstrap
resampling. However the Bayesian framework is more
suitable because it inherently allows to model uncertainty

and to incorporate a priori knowledge on the model pa-
rameters, as illustrated in Mélèse et al. (2018) for IDF
relationships.

In this article, a statistical model is proposed to assess
the severity of a storm at a given location and the related
uncertainties. We consider the IDAF model of Mélèse
et al. (2019) combined with a Gumbel distribution, in the
Bayesian framework. We estimate this model at two grid
points of a high-resolution reanalysis database covering a
mountainous region of southern France. The two con-
sidered locations lie a few kilometers apart—one lying
over the mountain crest while the other lies in the
mountain plateau. This proximity allows us to study the
spatial variability of the severity at short distance on ei-
ther side of the crest.

The data are presented in section 2. The Bayesian
framework of extreme value IDAF relationships al-
lowing us to assess the storm severity and the related
uncertainties is developed in section 3. Results are pre-
sented in section 4, with a particular focus on the severity
of the November 2011 event at the two locations.

## 2. Data

The region of study is centered on the Massif Central
range, a mountainous area of southern France that is
under Mediterranean climatic influence. We focus our
analysis on two locations: one is located at the Mont
Aigoual peak (1565 m MSL), in the southern part of the
Massif Central crest. The other one lies in the Massif
Central plateau, near the city of Meyrueis (800 m MSL),
15 km northwest of Mont Aigoual (see the white crosses in
Fig. 1). Although close in distance, these two locations
feature contrasting hydrometeorology. The Massif Central
crest delineates a topographical barrier protecting the
Massif Central plateau from Mediterranean moist winds
bringing extreme rainfall, giving usually much larger rain-
fall amounts over the crest than in the plateau (Mélèse
et al. 2018, 2019; Blanchet et al. 2016; Ceresetti et al. 2010).

We use high-resolution rainfall reanalysis data based on
radar–rain gauge merging, as described in Boudevillain
et al. (2016). These reanalyses are available at 1-km² spatial
resolution and at 1-h temporal resolution for 50 significant
rain events for the autumns 2008–15. The duration of these
events varies between a few hours to more than 8 successive
days. By studying extremes based on only eight years of
data, we are aware that interpretation of the results appeals
for vigilance. In particular, the large return period estima-
tions (say larger than 20 years) need to be considered
with prudence due to the large extrapolation they induce.
However the contribution of this article is primary a “proof
of concept” of a statistical framework allowing to estimate
the multiscale severity of a storm and the related
<figure><img src="image_2.png" /><figcaption>FIG. 1. Accumulated rainfall depth (mm km⁻¹) from 0600 UTC 1 Nov to 0000 UTC 6 Nov 2011 in the study region. The two white crosses indicate the localization of the two considered pixels, Meyrueis and Mont Aigoual. The white triangles show the main peaks of the Massif Central crest, and the white squares show the main towns.</figcaption></figure>

uncertainties. Although we use the best high-resolution spatiotemporal data available in our region, more reliable results would obviously benefit from longer time series.

Considering the 46 durations 3, 4, ..., 48 h and the 23 areas 1, 3², ..., 45² km², we compute the time series of areal rainfall centered on each of the two pixels, for each pair of duration and area, considering squared surfaces centered on the focus pixels for the areal aggregation. Then two databases are constructed for each pixel. DB1 contains the annual maximum intensities for each pair of duration and area. This database will be used for the estimation of the IDAF model (section 3). For computational reasons and owing to the fact that very dependent scales bring almost the same information, we restrict this database to 9 durations (3, 4, 6, 8, 12, 16, 24, 36, and 48 h) and 10 areas (1, 9, 25, 81, 169, 361, 625, 1089, 1521, and 2025 km²), giving for each pixel and each year a vector of 90 annual maxima. DB2 contains the event maximum intensities for all spatiotemporal scales (without restriction), giving for each pixel and each event a vector of 1058 event maxima. This second database will be used to compute the multiscale severity of the 50 rain events, based on the estimated IDAF model (section 4). A particular focus will be made on the event that occurred during 1–6 November 2011, which resulted in the largest rainfall accumulation over the Massif Central crest during the observation period.

### 3. Multiscale severity modeling

The model of IDAF relationships considered in this article has been developed in De Michele et al. (2001) in a frequentist framework and used in Panthou et al. (2014). It is a particular case of the model developed in Mélèse et al. (2019) for the study region, using the same data. However, a Bayesian version is considered here, in which the IDAF parameters are considered as random variables rather than scalars. The proposed model is also an extension of the Bayesian IDF model considered in Mélèse et al. (2018) and Van de Vyver (2015) when adding the space component. For the sake of clarity, in the following, we write random variables with bold symbols to distinguish them from scalars. We first present the model in the frequentist framework, as in Mélèse et al. (2019), and then we develop the Bayesian case and its inference.

#### a. Extreme value IDAF relationships: Frequentist framework

The considered model applies to the annual maxima of DB1. The frequentist form of the model was deeply studied in Mélèse et al. (2019). Below we give its key elements. The model relies on two assumptions. The first one follows extreme value theory (Coles 2001). It is used to model the cumulative distribution function (CDF) of annual maximum intensity for some reference area $A_0$ and some reference duration $D_0$, taken in this article equal to the finest scale of the database, $A_0 = 1$ km² and $D_0 = 3$ h. The founding theorem of extreme value theory (Coles 2001, chapter 3) states that if independent and identically distributed data are blocked into sequences of observation and if each block is long enough, then the CDF of block maxima is approximately the generalized extreme value (GEV) distribution. The Gumbel distribution is a particular case of the GEV distribution in which the so-called shape parameter is fixed to 0. Although less flexible in general, section 4a will show that the Gumbel distribution is preferable for our data owing to the difficulty in estimating the GEV shape parameter based on short records (here 8 years). Thus, we suppose that the random variable of annual maximum intensity for the duration and area of reference $M(D_0, A_0)$ follows a Gumbel distribution, that is,

$$ \mathrm{Pr}\{M(D_0, A_0) < x\} = \exp\left[-\exp\left(-\frac{x - \mu_0}{\sigma_0}\right)\right], \quad (1) $$

where $\mu_0$ is a scalar representing the location parameter of the Gumbel distribution and $\sigma_0$ is the nonnegative scalar representing the scale parameter of the Gumbel distribution.
The second hypothesis relies on space-time structural properties of the rain fields. Assuming space-time scaling relationships, Mélèse et al. (2019) derived the CDF of annual maximum rainfall intensity for any pair of duration and area $M(D, A)$ from the CDF of annual maximum rainfall intensity for the duration and area of reference $M(D_0, A_0)$, as

$$ \Pr\{M(D, A) < x\} = \Pr\left\{\left(\frac{D}{D_0}\right)^{-H} \frac{s(D, A)}{s(D_0, A_0)} M(D_0, A_0) < x\right\}, \quad (2) $$

where $0 < H < 1$ is the temporal scaling exponent and $s(D, A) > 0$ models the spatiotemporal scaling. The model of De Michele et al. (2001) considers, up to a change of notation, that $s(D, A) = (1 + \omega D^{-\beta}A^{\alpha})^{\gamma}$, where $\omega, \beta, \alpha$, and $\gamma$ are scalars. However, Mélèse et al. (2019) showed that this is overparametrized for our data since there is no loss of generality in considering $\gamma = 1$ due to compensation effects between $\omega, \beta, \alpha$, and $\gamma$. The authors also showed that a more flexible model is to consider $s(D, A) = 1 + (\omega D^{-\beta} + \omega_2 D^{-\beta_2})A^{\alpha}$ in order to model that different scaling relationships may govern areal rainfall depending on the duration. However, as will be seen in section 4a, the additional term $\omega_2 D^{-\beta_2}$ reveals to be unnecessary for the two pixels of this study. Thus, we consider in this study

$$ s(D, A) = 1 + \omega D^{-\beta} A^{\alpha}. \quad (3) $$

Combining Eqs. (1)–(3) implies that the random variable of annual maximum rainfall intensity for any duration $D$ and any area $A$ follows a Gumbel distribution, that is,

$$ \Pr\{M(D, A) < x\} = \exp\left\{-\exp\left[-\left(\frac{D}{D_0}\right)^H \left(\frac{1 + \omega D^{-\beta} A^{\alpha}}{1 + \omega D_0^{-\beta} A_0^{\alpha}}\right)^{-1} \frac{x}{\sigma_0} + \frac{\mu_0}{\sigma_0}\right]\right\}, \quad (4) $$

where $\theta = (\mu_0, \sigma_0, H, \omega, \beta, \alpha)$ is the set of model parameters, with $H \ge 0, \sigma_0 > 0$ and $1 + \omega D^{-\beta}A^{\alpha} > 0$ for all $D, A$.

Let consider an event maximum $i$ in DB2 associated with a duration $D$ and an area $A$. Its return period is obtained from Eq. (4) combined with the expression of Stedinger et al. (1993) linking the return period of large values to the return period of annual maxima, giving

$$ T_{D,A,i} = \left\{ \exp\left[-\left(\frac{D}{D_0}\right)^H \left(\frac{1 + \omega D^{-\beta} A^{\alpha}}{1 + \omega D_0^{-\beta} A_0^{\alpha}}\right)^{-1} \frac{i}{\sigma_0} + \frac{\mu_0}{\sigma_0}\right] \right\}^{-1}. \quad (5) $$

Parameter $T_{D,A,i}$ can be interpreted as the mean number of years between two consecutive rainfall intensities at scales ($D, A$) that exceed the intensity $i$. Following Ramos et al. (2005), we define the multiscale severity of a given event as the set of return periods associated with all event maxima, when considering the continuum of durations and areas. Note that a slightly different definition is considered in Cattiaux and Ribes (2018), where the minimum probability of exceedance is considered rather than the return period of the event maximum.

b. Extreme value IDAF relationships: Bayesian framework

The Bayesian framework differs from the frequentist framework in that the model parameters are considered as random variables rather than as scalars. As will be seen below, Bayesian inference gives a direct and flexible estimation of uncertainty through probability density functions. Mélèse et al. (2018) illustrated this flexibility with respect to the frequentist case for IDF relationships. The above Eqs. (1)–(4) still hold in the Bayesian case but the probabilities on the left-hand sides are conditional to the random variable of model parameters $\theta$ being equal to some value $\theta = (\mu_0, \sigma_0, H, \omega, \beta, \alpha)$. The return period is also given by Eq. (5) but the equality applies to random variables.

For the sake of conciseness, we write $\underline{m}$ the set of annual maximum rainfall intensities for all durations and all areas, that is, $\underline{m} = \{m_{D,A,i}; D \in \mathcal{D}, A \in \mathcal{A}, y \in \mathcal{Y}\}$, where $\mathcal{D}$ is the set of durations, $\mathcal{A}$ is the set of areas, and $\mathcal{Y}$ is the set of years. In the Bayesian framework (Gelman et al. 2004), inference requires estimating the density of the parameters knowing the data, that is, $f(\theta|\underline{m})$, called the posterior density. The Bayes formula states that

$$ f(\theta|\underline{m}) = \frac{f(\underline{m}|\theta)f_{\theta}(\theta)}{\int_{\theta} f(\underline{m}|\theta)f_{\theta}(\theta)d\theta}, \quad (6) $$

where $f(\underline{m}|\theta)$ is the likelihood function and $f_{\theta}$ is the density of $\theta$ called the prior density. We suppose independence of the model parameters, that is,

$$ f_{\theta}(\theta) = f_{\mu_0}(\mu_0) f_{\sigma_0}(\sigma_0) f_H(H) f_{\omega}(\omega) f_{\beta}(\beta) f_{\alpha}(\alpha), \quad (7) $$

where $f$ denote univariate densities. In this study, the six chosen prior densities are weakly informative, using a uniform density for $H$ and Gaussian densities for $\mu_0, \sigma_0, \omega, \beta, \alpha$ (see details in section a of the appendix). Under the assumptions that: (i) annual maxima are independent from one year to another and (ii) annual maxima
of a given year at different durations and areas are independent, the expression of the likelihood function $f(\underline{m}|\underline{\theta})$ is given by

$$f_{\text{ind}}(\underline{m}|\underline{\theta}) = \prod_{y \in \mathcal{Y}} f_y(\underline{m}|\underline{\theta}) = \prod_{y \in \mathcal{Y}} \prod_{A \in \mathcal{A}} \prod_{D \in \mathcal{D}} f_{D,A}(m_{D,A,y}|\underline{\theta}), \quad (8)$$

where $f_y$ is the density associated to all the maxima of year $y$ and $f_{D,A}$ is the Gumbel density associated with Eq. (4). However the assumption ii) above is likely to be misspecified. As a consequence, the likelihood $f_{\text{ind}}$ is likely not to be the true probability of the observed data, and the posterior density in Eq. (6) under Eq. (8) can be misleading. Two adjustments of composite likelihoods such as the independence likelihood $f_{\text{ind}}$ have been proposed in Pauli et al. (2011) and Chandler and Bate (2007) in order to retrieve some of the desirable properties of the full likelihood: the so-called *magnitude adjustment* and the *curvature adjustment*. Ribatet et al. (2012) considered both of them for extremes in a Bayesian framework. Van de Vyver (2015) compared them for IDF relationships and found few differences. A separate analysis reveals that for our data the curvature adjustment induces a slight bias in the quantile estimates of Mont Aigoual, possibly due to the fact that $1 + \omega D^{-\beta}A^{\alpha}$ in Eq. (4) can be negative for some $\theta$ even in the close neighborhood of $\hat{\theta}_{\text{ind}}$. Therefore only the magnitude adjustment is considered hereinafter. It replaces the independence likelihood $f_{\text{ind}}$ by the adjusted likelihood:

$$f_{\text{adj}}(\underline{m}|\underline{\theta}) = \{f_{\text{ind}}(\underline{m}|\underline{\theta})\}^k, \quad (9)$$

where $k$ is a scalar whose computation is explained below. Plugging Eq. (9) into Eq. (6) gives the adjusted posterior density

$$f_{\text{adj}}(\underline{\theta}|\underline{m}) \propto f_{\text{adj}}(\underline{m}|\underline{\theta})f_{\theta}(\underline{\theta}), \quad (10)$$

where “$\alpha$” means “proportional to.”

Markov chain Monte Carlo (MCMC) techniques provide a way of simulating from complex distribution with untrackable or unknown normalizing constant—such as the adjusted posterior density $f_{\text{adj}}(\underline{\theta}|\underline{m})$ in Eq. (10)—by simulating Markov chains that have the target distributions as their stationary distributions. The posterior density can be obtained from the simulated samples at convergence of the Markov chains. One of the most popular MCMC techniques is the Gibbs sampling algorithm. Consider the vector of parameters $\theta = (\mu_0, \sigma_0, H, \omega, \beta, \alpha)$ that for ease and generality of presentation we relabel $\theta = (\theta_1, \dots, \theta_n)$ with $n = 6$. At each step of the Gibbs sampler, a new sample of $\theta_j$ is obtained iteratively, cycling through the $n$ parameters

from $j = 1$ to $j = n$, by drawing from the full adjusted conditional density

$$f_{\text{adj}}(\underline{\theta}_j|\tilde{\underline{\theta}}_{-j}, \underline{m}) \propto \{f_{\text{ind}}(\underline{m}|\underline{\theta}_j, \tilde{\underline{\theta}}_{-j})\}^k f_{\theta_j}(\underline{\theta}_j), \quad (11)$$

where $\tilde{\underline{\theta}}_{-j}$ is the current sample for the other parameters. As in our case the full conditional densities $f_{\text{adj}}(\underline{\theta}_j|\tilde{\underline{\theta}}_{-j}, \underline{m})$ cannot be sampled directly, a Metropolis-Hasting step is substituted. The details of the algorithm are given in section b of the appendix. Two cases are considered for the scalar $k$ involved in the likelihood adjustment (Ribatet et al. 2012; Van de Vyver 2015). In the overall adjustment, $k$ is computed before running the Gibbs algorithm as

$$k = n/\text{tr}[\mathbf{I}(\hat{\theta}_{\text{ind}})^{-1} \mathbf{V}(\hat{\theta}_{\text{ind}})], \quad (12)$$

with $\hat{\theta}_{\text{ind}}$ the set of parameters maximizing the independence likelihood (8), tr the trace of the matrix and $\mathbf{I}(\hat{\theta}_{\text{ind}})$ and $\mathbf{V}(\hat{\theta}_{\text{ind}})$ the $n \times n$ matrices $\mathbf{I}(\underline{\theta}) = E[-\nabla^2 \log f_{\text{ind}}(\underline{m}|\underline{\theta})]$ and $\mathbf{V}(\underline{\theta}) = \text{cov}[\nabla \log f_{\text{ind}}(\underline{m}|\underline{\theta})]$ evaluated in $\hat{\theta}_{\text{ind}}$, which are approximated by

$$\mathbf{I}(\underline{\theta}) = - \sum_{y \in \mathcal{Y}} \frac{\partial^2 \log f_y(\underline{m}|\underline{\theta})}{\partial \underline{\theta} \partial \underline{\theta}^T}, \quad (13)$$

$$\mathbf{V}(\underline{\theta}) = \sum_{y \in \mathcal{Y}} \frac{\partial \log f_y(\underline{m}|\underline{\theta}) \partial \log f_y(\underline{m}|\underline{\theta})}{\partial \underline{\theta} \partial \underline{\theta}^T}. \quad (14)$$

In the *adaptative adjustment*, $k$ is updated at each step of the Gibbs sampler and for each parameter $\theta_j$ as

$$k = I_j(\hat{\theta}_{\text{ind},j})/V_j(\hat{\theta}_{\text{ind},j}), \quad (15)$$

where $\hat{\theta}_{\text{ind},j}$ is the scalar maximizing the independence likelihood in Eq. (8) with respect to $\theta_j$, with the other parameters holding fixed to their current values $\hat{\theta}_{-j}$, and $I_j(\hat{\theta}_{\text{ind},j})$ and $V_j(\hat{\theta}_{\text{ind},j})$ are the two scalars

$$I_j(\underline{\theta}) = - \sum_{y \in \mathcal{Y}} \frac{\partial^2 \log f_y(\underline{m}|\underline{\theta}_j, \tilde{\underline{\theta}}_{-j})}{\partial \underline{\theta}_j^2}, \quad (16)$$

$$V_j(\underline{\theta}) = \sum_{y \in \mathcal{Y}} \left[ \frac{\partial \log f_y(\underline{m}|\underline{\theta}_j, \tilde{\underline{\theta}}_{-j})}{\partial \underline{\theta}_j} \right]^2, \quad (17)$$

evaluated in $\hat{\theta}_{\text{ind},j}$. These two adjustments will be compared in section 4.

Given the MCMC samples, samples of the posterior density $f(T_{D,A,i}|\underline{m})$ are obtained by applying Eq. (5) to the samples of $\mu_0, \sigma_0, H, \omega, \beta, \alpha$. This gives an estimate of the return period posterior density $f(T_{D,A,i}|\underline{m})$. When needed, the latter will be summarized by its mode,
<table><caption>TABLE 1. Composite likelihood information criteria for Mont Aigoual and Meyrueis.</caption><thead><tr><th>s(D, A)</th><th>EV distribution</th><th>Mont Aigoual</th><th>Meyrueis</th></tr></thead><tbody><tr><td>(1 + ωD<sup>-β</sup>A<sup>α</sup>)<sup>γ</sup></td><td>Gumbel</td><td>3571.3</td><td>2684.3</td></tr><tr><td>(1 + ωD<sup>-β</sup>A<sup>α</sup>)<sup>γ</sup></td><td>GEV</td><td>3596.0</td><td>—<sup>a</sup></td></tr><tr><td>1 + (ωD<sup>-β</sup> + ω<sub>2</sub>D<sup>-β<sub>2</sub></sup>)A<sup>α</sup></td><td>Gumbel</td><td>3570.1</td><td>—<sup>b</sup></td></tr><tr><td>1 + (ωD<sup>-β</sup> + ω<sub>2</sub>D<sup>-β<sub>2</sub></sup>)A<sup>α</sup></td><td>GEV</td><td>3596.7</td><td>—<sup>b</sup></td></tr><tr><td>1 + ωD<sup>-β</sup>A<sup>α</sup></td><td>Gumbel</td><td>3569.7</td><td>2681.3</td></tr><tr><td>1 + ωD<sup>-β</sup>A<sup>α</sup></td><td>GEV</td><td>3596.3</td><td>—<sup>a</sup></td></tr></tbody></table>

which is interpreted as the most likely mean number of years between two consecutive exceedances of the level $i$, together with the 2.5% and 97.5% quantiles, respectively noted $T_{D,A,i}^{\text{mode}}$, $T_{D,A,i}^{\text{low}}$, and $T_{D,A,i}^{\text{up}}$. Likewise, the multiscale severity of an event will be summarized by the set of ($T_{D,A,i}^{\text{mode}}$, $T_{D,A,i}^{\text{low}}$, $T_{D,A,i}^{\text{up}}$) when $D$ and $A$ range the continuum of scales.

# 4. Results

## a. Model selection (frequentist framework)

We temporarily consider the three different models for $s(D, A)$ of section 4a and evaluate which one is the most suitable for our data, together with the choice of the extreme value distribution (either Gumbel or GEV), in the frequentist framework. We consider the independence likelihood of Eq. (8) and compute the most likely set of parameters $\hat{\theta}_{ind}$ in each case. Model selection is based on comparing information criteria. However since the independence likelihood is not the “true” likelihood of the data (due to dependence among areas and durations), the regular information criteria such as Akaike’s information criterion (AIC) cannot be used. According to Varin and Vidoni (2005), the analog of the AIC for composite likelihoods (among which the independence likelihood) is the composite likelihood information criterion (CLIC). It is given by (Varin and Vidoni 2005; Varin 2008)

$$ \text{CLIC} = -2 \log f_{\text{ind}}(m|\hat{\theta}_{\text{ind}}) + 2\text{tr}[\mathbf{I}(\hat{\theta}_{\text{ind}})^{-1}\mathbf{V}(\hat{\theta}_{\text{ind}})]. \quad (18) $$

For two competing models, lower values of CLIC correspond to better quality models. Table 1 confirms that the model $s(D, A) = 1 + ωD^{-β}A^{α}$ with the Gumbel distribution—as assumed in the article—is preferable for Mont Aigoual. Note that the GEV shape parameter for Mont Aigoual is estimated to -0.05 for the three models, which is very close to 0 corresponding to a Gumbel distribution. For Meyrueis, numerical issues

prevent from computing the CLIC apart for the models $s(D, A) = (1 + ωD^{-β}A^{α})^{γ}$ and $s(D, A) = 1 + ωD^{-β}A^{α}$ in the Gumbel case, with a preference for the later according to CLIC. The GEV shape parameter for Meyrueis is estimated to -0.30, corresponding to a bounded (reversed Weibull) distribution but the $\mathbf{I}(\hat{\theta}_{ind})$ and $\mathbf{V}(\hat{\theta}_{ind})$ matrices cannot be computed (see the caption of Table 1). Furthermore, despite the shape parameter quite different from 0, the estimated quantiles are actually very similar in the GEV and Gumbel case (absolute difference lower than 1 mm h<sup>-1</sup>; not shown). Based on this observation and given that the $\mathbf{I}(\hat{\theta}_{ind})$ and $\mathbf{V}(\hat{\theta}_{ind})$ matrices are necessary to adjust the likelihood in the Gibbs sampler, we select for the rest of this study the model $s(D, A) = 1 + ωD^{β}A^{α}$ with the Gumbel distribution, for both pixels.

## b. MCMC monitoring

The IDAF model developed in section 3 is estimated at the two considered locations based on the annual maxima of DB1. For each pixel, we run the Gibbs algorithm four times, giving four MCMC chains whose convergence is monitored using the $\hat{R}$ convergence criteria of Gelman et al. (2014, chapter 6). We consider that convergence is reached for $\hat{R} < 1.005$, which is obtained after 200 000 iterations. The burn-in period is set to the first half of iterations. Every 100th iteration of the remaining 100 000 iterations of the four MCMC chains is considered for estimating the posterior densities in order to reduce the dependence within the sample. The posterior density estimates are obtained from these 4000 iterations.

## c. Model validation

We assess goodness-of-fit of the Bayesian Gumbel-IDAF model by comparing the empirical and modeled quantiles. The empirical quantiles at a given scale ($D, A$) are the sorted annual maxima. Let $m_{D,A,(r)}$ be the $r$th smallest maxima in the annual maxima series for duration $D$ and area $A$ at a given pixel, for $r = 1, \dots, 8$.

<footnote><sup>a</sup> CLIC cannot be computed because $f_{ind}$ is undefined ($1 + ωD^{-β}A^{α} < 0$) in the neighborhood of $\hat{\theta}_{ind}$.</footnote>

<footnote><sup>b</sup> CLIC cannot be computed because $β_2$ is very large (>200), giving 0 values in the $\mathbf{I}(\hat{\theta}_{ind})$ matrix. Note that when $β_2$ is very large, $D^{-β_2} \approx 0$, giving $s(D, A) \approx 1 + ωD^{-β}A^{α}$.</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 2. Q–Q plot (mm h⁻¹) for (top) Meyrueis and (bottom) Mont Aigoual for all spatiotemporal scales (D, A). The blue line shows the posterior modes of the quantile densities. The red lines show the 95% credibility intervals.</figcaption></figure>

This quantile is associated with a nonexceedance probability $p_r$. Different definitions are possible for $p_r$; here we use Blom’s plotting positions (Blom 1958) $p_r = (r - 0.375)/(8 + 0.25)$, giving a nonexceedance probability of 0.076 for the smallest maxima and 0.924 for the largest maxima. The estimated quantile associated with this probability is obtained by inverting $Pr\{M(D, A)|\theta\} = p_r$ using Eq. (4). Plugging the MCMC samples of the parameters $\theta$ in Eq. (4) gives a sample of posterior quantiles associated with $p_r$, and hence an estimate of its posterior density. Figure 2 compares the empirical quantiles with the posterior quantiles, for all possible pairs of scales, summarizing the quantile posterior density using the posterior mode together with the 95% credibility interval. Three cases are considered: without adjustment of the likelihood, or with either the overall or adaptative magnitude adjustment. Figure 2 shows that, whatever case, the posterior modes give globally unbiased estimates of the empirical quantiles. In particular, there is no systematic over or underestimation of the quantiles despite the use of a Gumbel

rather than a GEV distribution. The three cases give very similar posterior modes for Meyrueis, however we note a larger variability of the posterior modes for Mont Aigoual with the adaptative adjustment. Even more striking is the difference in the credibility intervals when adjusting or not the likelihood. Not adjusting the likelihood clearly underestimate the uncertainty, as also shown in Van de Vyver (2015) for IDF relationships. Fifty-eight percent of the empirical quantiles for Meyrueis fall outside the 95% credibility intervals when the likelihood is not adjusted and 63% for Mont Aigoual. For the two pixels, the adaptative adjustment gives slightly larger credibility intervals than the overall adjustment, increasing mainly the upper bound of the largest quantiles. Of the empirical quantiles for Meyrueis and Mont Aigoual, 14% and 24%, respectively, fall outside the 95% credibility intervals with the overall adjustment, versus 5% and 10% with the adaptative adjustment (recall that under the true model, 5% of the quantiles are expected to fall outside the 95% credibility interval). This leads us to select
<figure><img src="image_4.png" /><figcaption>FIG. 3. Posterior density estimates of the 2011 event return periods in Mont Aigoual and Meyrueis at (left) *D* = 3 h, *A* = 81 km² and (right) *D* = 24 h, *A* = 1089 km². The dashed vertical lines show the mode of the posterior density. The colored segments on the *x* axis show the lower and upper bounds of the 95% credibility intervals.</figcaption></figure>

the adaptative adjustment for the rest of this study. Although Van de Vyver (2015) found almost no difference between the two adjustments for IDF relationships, the fact that the adaptative adjustment gives better uncertainty estimates is not surprising since it is more in line with the Gibbs sampling “spirit,” which is to sample the one dimensional density of $\theta_j$ given the current value of $\tilde{\theta}_{-j}$ (Ribatet et al. 2012).

### d. Multiscale severity of the 2011 event

A particular focus is made on the event of DB2 that gave the largest rainfall accumulation over the Massif Central crest. This event occurred during 1–6 November 2011. Although several deep convective systems were triggered during 1, 4, and 5 November in the foothills, the great majority of this event corresponds to quasi-continuous orographic rainfall in time and space over the Massif Central crest from 2 to 6 November (Boudevillain et al. 2016). Figure 1 displays the map of the accumulated rainfall depth during this event. More than 700 mm rainfall was recorded at the South of Serre Croix de Bauzon and Mont Aigoual. This event is a typical case of shallow convection characterized by north–south rainbands triggered by the interaction of moist southern winds and details of the topography (Miniscloux et al. 2001; Anquetin et al. 2003; Godart et al. 2009).

Figure 3 compares the return period posterior density estimates of this event. For illustration purposes, we focus here on two scales that are relevant from a hydrological point of view (Marchi et al. 2010), considering (3 h, 81 km²) and (24 h, 1089 km²). For either pixel, the

posterior density varies greatly when passing from one spatiotemporal scale to another, which was obviously expected since rainfall neither accumulates uniformly in time nor in space. For both scales, Mont Aigoual and its neighborhood received twice as much rainfall as Meyrueis, illustrating the mask effect induced by the crest that protects the Massif Central plateau from the moist southern winds bringing rainfall extremes. Despite this factor of 2, the return period posterior densities are quite similar at both locations. The 2011 event featured long lasting orographic rainfalls leading to quite ordinary rainfall intensities for short durations (return periods of a few years) but much more extraordinary accumulations at the daily scale (return periods of several tens of years).

Figure 3 also reveals that whatever scale and pixel, the posterior densities are asymmetric with heavier right tails than the left tails, even more in Mont Aigoual than in Meyrueis. To quantify this, we consider the ratio of asymmetry defined as $e = (T_{D,A,i}^{\text{up}} - T_{D,A,i}^{\text{mode}})/(T_{D,A,i}^{\text{mode}} - T_{D,A,i}^{\text{low}})$. Values of $e = 1, e > 1,$ and $e < 1$ respectively correspond to symmetric, right-skewed, and left-skewed posterior densities. The asymmetry ratio $e$ is 8 (6) in Mont Aigoual (Meyrueis) for the (3 h, 81 km²) scale and 16 (9) for the (24 h, 1089 km²) scale. This implies that, for both scales and both locations, there is actually a much greater likelihood for the true 2011 return periods to be larger than $T^{\text{mode}}$ rather than lower. Estimating the return period as its most likely values $T^{\text{mode}}$, omitting uncertainties, is likely to underestimate the true return period.

As a way of generalizing the previous results for the continuum of scales, we visualize in Fig. 4 the multiscale
<figure><img src="image_5.png" /><figcaption>FIG. 4. Severity of the 2011 event in Meyrueis and Mont Aigoual: (a) event maximum, (b) posterior mode of the return period densities, (c) lower bound of the 95% credibility interval, and (d) upper bound of the 95% credibility interval. Note that the color scale changes from one plot to another.</figcaption></figure>
maxima together with the multiscale severity of the 2011 event, in the vein of the severity diagrams of Ramos et al. (2005) but considered here at a given location rather than at regional scale. The maximum intensity diagrams of the two locations confirm the mask effect induced by the topography, with decreasing maximum rainfall intensity as the area increases in Mont Aigoual, whereas the contrary applies to Meyrueis.

Keeping in mind the shortness of the database, we comment below on the main patterns of the severity diagrams, interpreting the largest return period estimates in a rough way rather than commenting on the exact numbers. Comparing the diagrams of maximum intensity and severity mode shows that the spatiotemporal scales featuring the largest maxima do not correspond to those featuring the largest return periods. In Mont Aigoual, the maximum rainfall intensities produced during the 2011 event were quite ordinary for durations of 3–4 h, with most likely return periods around a few years. On the contrary, they were very exceptional for durations ranging 16–36 h over areas extending from 10 to 250 km², with most likely return periods of several tens of years. This shows that, at a given location, a storm can be both ordinary and exceptional depending on the considered spatiotemporal scales, as also shown in Ramos et al. (2005). The most severe spatiotemporal scales (i.e., corresponding to the largest return periods) for the 2011 event in Meyrueis are 10–30 h and 1000–2000 km². In Mont Aigoual, similar temporal scales are found (16–30 h) but they are associated with much lower areas (10–200 km²) due to the topographical effect intensifying rainfall over the crest.

As expected, the patterns of the lower and upper bounds of the credibility intervals are broadly similar to that of the mode since larger return periods tend to come with larger uncertainties, despite some nonlinearities. For both locations, there is a larger uncertainty on the right side across all scales but particularly for the most severe scales and for Mont Aigoual where the asymmetry ratio is up to 63 (up to 37 in Meyrueis).

### e. Generalization to other events

The asymmetry of the return period densities shown above for the 2011 event is actually generalizable to all the significant event maxima across all scales. To show this, we consider the event maxima of DB2 whose posterior modes are larger than 2 years, excluding the 2011 event that was already described in the previous section. There are 1270 such cases for Mont Aigoual and 1895 for Meyrueis. Figure 5 shows the asymmetry ratio *e* for each of these event maxima. There is no left-skewed return period density (i.e., all asymmetry ratios are

larger than 1). The 15 event maxima of Mont Aigoual whose mode exceeds 10 years show an average *e* of 11.

Figure 6 shows the posterior modes of the severity diagrams of the four secondmost severe events in Mont Aigoual. The 2010 and 2008 events (Figs. 6a,b) are cases of shallow convection organized in north-south rainbands, as in the 2011 event (see Fig. 1). The 2015 and 2014 events (Figs. 6c,d) show situations of deep convection during which the circulation of the mesoscale convective system (MCS) is modified by the interaction with the topography at different scales. This interaction more or less stabilizes the MCS upstream from the crest, over the foothills (Ducrocq et al. 2008). Considering the pattern of the severity diagrams (i.e., the colors) rather than the absolute values of the return periods, Fig. 6 reveals the variety of severity patterns both within a given class of precipitation (deep convection or shallow convection) and between the two classes. Although the pattern of the severity diagram of the 2010 event may appear somewhat similar to that of the 2011 event (see Fig. 4), it actually shows a case of severe precipitation over much shorter durations. Broadly speaking, in both Meyrueis and Mont Aigoual, the spatiotemporal “dynamics” of the 2010 event (as represented by the return period) was similar to that of the 2011 event but during a much shorter time. The 2008 event, although corresponding to the same class of shallow banded convection, shows a much more different pattern of severity. It was indeed quite ordinary in Mont Aigoual and its

<figure><img src="image_6.png" /><figcaption>FIG. 5. Boxplot of the ratio of asymmetry *e* for Mont Aigoual and Meyrueis for all event maxima with $T_{D,Aj}^{\text{mode}} \ge 2$ years, excluding the 2011 event. The upper and lower red points show the quantiles of order 0.975 and 0.025, respectively. The blue horizontal line at 1 corresponds to a symmetric density.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIG. 6. Accumulated rainfall depth (mm) and severity diagrams (posterior modes) of the four second-most severe events in Mont Aigoual: (a) from 0600 UTC 3 Oct to 0000 UTC 5 Oct 2010, (b) from 0600 UTC 31 Oct to 0600 UTC 3 Nov 2008, (c) from 0600 UTC 11 Sep to 0600 UTC 18 Sep 2015, and (d) from 0600 UTC 15 Sep to 0000 UTC 21 Sep 2014.</figcaption></figure>
neighborhood at subdaily to daily scale but it was more unusual at all spatial scales for 24-48-h durations. Rainfall intensities were basically less intense but more spread in time and space around Mont Aigoual than during the 2010 and 2011 events. The two cases of deep convection show also very different severities in both Mont Aigoual and Meyrueis. In Mont Aigoual, the largest severities of the 2015 event are found for quite small areas (10-300 km²) and medium durations (6-24 h), roughly like the 2011 event of Fig. 4. However a noteworthy difference is that the medium severities (relatively to the largest values) extend down to 3 h and up to 2000 km² due to the stabilization of the MCS upstream from the crest. On the contrary, the 2015 event shows almost uniform severities over all spatiotemporal scales. In Meyrueis, the daily scale and large areas feature among the largest severities, which is again the signature of the MCS stabilization upstream from the crest. Although not exhaustive, these examples illustrate that the multiscale severity diagram is a powerful tool to visualize the “specialness” of a given event at a given location. Neither an accumulation map nor a couple of return period estimates can provide such a comprehensive description of how “special” an event was across the scales.

**5. Conclusions**

In this article we assessed the multiscale severity of a given storm at a given location. A Bayesian framework was proposed to compute the severity and the associated uncertainties at multiple scales. The multiscale severity of 50 events were computed at two nearby locations of the Massif Central, using a rain gauge-radar reanalysis database. The severity diagram revealed to be a very useful tool to efficiently visualize the severity of a storm at a given location or to compare the severities at different locations and/or of different storms. Three main results were found. First, the study revealed right-skewed posterior densities of return periods across all scales, and particularly for the largest return periods, evidencing the benefit of considering uncertainty in order to avoid a systematic underestimation of the frequency of occurrence of the severe events. Second, we showed the variability of the storm return period at a given location depending on the spatiotemporal scales, showing the benefit of assessing and visualizing the severity over the continuum of durations and areas rather than at predetermined spatiotemporal scales. Finally, we showed the strong spatial variability of the severity at short distance. This variability is partly due to the inherent rainfall variability but it is exacerbated here by the topography of the region. This result emphasizes the

merit of using high-resolution radar data rather than a sparser rain gauge network.

Despite a significant progress brought in the area of rainfall severity assessment, a limitation of the present study is that the considered spatial scales of aggregation, which are squared, do not coincide with the hydrological risk that depends on the shape of watershed over which rainfall accumulates. An improvement of this study will be to consider the severity over areas matching the watersheds. However a difficulty is that the watersheds of interest in the region extend from the Massif Central crest to the foothills. The latter undergoes deep convective systems whose positioning is unpredictable, making their sampling difficult within a short observation period (Mélèse et al. 2019). Methods based on stochastic storm transposition (Foufoula-Georgiou 1989; Wright et al. 2013; Zhou et al. 2019) reshuffling the reanalysis data to construct longer series of rain fields in the foothills seem an interesting idea to pursue.

*Acknowledgments.* The authors warmly thank Brice Boudevillain and Guy Delrieu for providing us the radar reanalysis data and for fruitful discussions on the hydro-meteorology of the region. The data used in this paper are maintained by OHMCV (https://ohmcv.osug.fr/). They can be downloaded at http://doi.osug.fr/public/OHMCV/OHMCV.REA.CEV.07-14.1.html.

# APPENDIX

## Technical Details on the MCMC Algorithm

### a. *Choice of the priors*

For generality purpose, we use weakly informative priors for each of the six parameters, using either Gaussian or uniform densities. The parameters of these densities are set at region scale, based on the classification of the areal rainfall structure obtained in Mélèse et al. (2019) as follows. First, we consider the frequentist version of the IDAF model defined by Eq. (4), which we estimate by maximum likelihood at 2149 pixels spread over the Massif Central and the foothills, including Mont Aigoual and Meyrueis pixels. Then, as a way of grouping together the pixels featuring similar extreme areal rainfall structure, we classify with *K* means the normalized curves of areal reduction factor (ARF) derived from Eq. (2) as the multiplicative factor allowing to pass from the CDF of annual maxima for any (*D*, *A*₀) to that for any (*D*, *A*). Finally we consider the classes associated with Mont Aigoual and Meyrueis pixels, which are not the same.

For the scaling parameter *H*, we use a uniform prior density between 0 and 1 as in Van de Vyver (2015) and
Mélèse et al. (2018) in an IDF framework. For the other parameters we use Gaussian prior densities. The Gaussian mean parameters are set to the averages of the maximum likelihood estimates of $\mu_0$, $\sigma_0$, $\omega$, $\beta$, $\alpha$ within the corresponding class, that is, respectively to 16.8, 7.1, -0.04, 0.12, 0.4 for Mont Aigoual and 8.1, 3.1, -0.12, -0.003, 0.4 for Meyrueis. The Gaussian standard deviations are set to twice the standard deviation of the maximum likelihood estimates of $\mu_0$, $\sigma_0$, $\omega$, and $\alpha$ within the corresponding class and four times for $\beta$, that is, respectively to 8.4, 3.9, 0.36, 0.28, 0.44 for Mont Aigoual and 2.02, 1.02, 0.74, 0.91, 0.73 for Meyrueis.

## b. Gibbs sampling algorithm

The algorithm proceeds as follows (Van de Vyver 2015):

1) Draw a starting point $\theta^{(0)}$ for which $f[\theta^{(0)}|\underline{m}]$ is defined and nonzero.

2) At each step $t = 1, 2, \dots$,

2.1) Set $\tilde{\theta} = \theta^{(t-1)}$.

2.2) For each $j = 1, \dots, n$,

(i) Draw $\epsilon_j$ from a Gaussian distribution with mean 0 and variance $\sigma_j^2$ and set $\theta_j^* = \theta_j^{(t-1)} + \epsilon_j$.

(ii) Derive the acceptance probability $a_j$:

$$ a_j = \min \left\{ 1, \frac{f_{\text{adj}}(\theta_j^*|\tilde{\theta}_{-j}, \underline{m})}{f_{\text{adj}}(\theta_j^{(t-1)}|\tilde{\theta}_{-j}, \underline{m})} \right\} \\ = \min \left\{ 1, \frac{f_{\text{adj}}(\underline{m}|\theta_j^*, \tilde{\theta}_{-j})f(\theta_j^*)}{f_{\text{adj}}(\underline{m}|\theta_j^{(t-1)}, \tilde{\theta}_{-j})f(\theta_j^{(t-1)})} \right\}. \quad (A1) $$

(iii) Accept the candidate $\theta_j^*$, that is, set $\tilde{\theta}_j = \theta_j^*$, with probability $a_j$.

2.3) Set $\theta^{(t)} = \tilde{\theta}$.

In our case, $\theta^{(0)}$ is set to the maximum likelihood estimate for the pixel. The standard deviations $\sigma_j$ are chosen so that the acceptance rate (i.e., the proportion of candidates that are accepted) is between 30% and 50% since the optimal acceptance rate is approximately 0.44 for one dimensional target distributions (see Gelman et al. 2014, chapter 12).

## REFERENCES

Anquetin, S., F. Minsicloux, J.-D. Creutin, and S. Cosma, 2003: Numerical simulation of orographic rainbands. *J. Geophys. Res.*, **108**, 8386, https://doi.org/10.1029/2002JD001593.

Blanchet, J., D. Ceresetti, G. Molinié, and J.-D. Creutin, 2016: A regional GEV scale-invariant framework for intensity duration frequency analysis. *J. Hydrol.*, **540**, 82–95, https://doi.org/10.1016/j.jhydrol.2016.06.007.

Blom, G., 1958: Statistical Estimates and Transformed Beta-Variables. Wiley, 176 pp.

Boudevillain, B., G. Delrieu, A. Wijbrans, and A. Confoland, 2016: A high-resolution rainfall re-analysis based on radar-raingauge merging in the Cévennes-Vivarais region, France. *J. Hydrol.*, **541A**, 14–23, https://doi.org/10.1016/j.jhydrol.2016.03.058.

Bougadis, J., and K. Adamowski, 2006: Scaling model of a rainfall intensity-duration-frequency relationship. *Hydrol. Processes*, **20**, 3747–3757, https://doi.org/10.1002/hyp.6386.

Cattiaux, J., and A. Ribes, 2018: Defining single extreme weather events in a climate perspective. *Bull. Amer. Meteor. Soc.*, **99**, 1557–1568, https://doi.org/10.1175/BAMS-D-17-0281.1.

Ceresetti, D., G. Molinié, and J.-D. Creutin, 2010: Scaling properties of heavy rainfall at short duration: A regional analysis. *Water Resour. Res.*, **46**, W09531, https://doi.org/10.1029/2009WR008603.

—, S. Anquetin, G. Molinié, E. Leblois, and J.-D. Creutin, 2012: Multiscale evaluation of extreme rainfall event predictions using severity diagrams. *Wea. Forecasting*, **27**, 174–188, https://doi.org/10.1175/WAF-D-11-00003.1.

Chandler, R. E., and S. Bate, 2007: Inference for clustered data using the independence loglikelihood. *Biometrika*, **94**, 167–183, https://doi.org/10.1093/biomet/asm015.

Coles, S., 2001: *An Introduction to Statistical Modeling of Extreme Values*. Springer, 209 pp.

—, and L. Pericchi, 2003: Anticipating catastrophes through extreme value modelling. *J. Roy. Stat. Soc.*, **52C**, 405–416, https://doi.org/10.1111/1467-9876.00413.

De Michele, C., N. T. Kottegoda, and R. Rosso, 2001: The derivation of Areal Reduction Factor of storm rainfall from its scaling properties. *Water Resour. Res.*, **37**, 3247–3252, https://doi.org/10.1029/2001WR000346.

Ducrocq, V., O. Nuissier, D. Ricard, C. Lebeaupin, and T. Thouvenin, 2008: A numerical study of three catastrophic precipitating events over southern France. II: Mesoscale triggering and stationarity factors. *Quart. J. Roy. Meteor. Soc.*, **134**, 131–145, https://doi.org/10.1002/qj.1999.

Foufoula-Georgiou, E., 1989: A probabilistic storm transposition approach for estimating exceedance probabilities of extreme precipitation depths. *Water Resour. Res.*, **25**, 799–815, https://doi.org/10.1029/WR025i005p00799.

Gelman, A., J. B. Carlin, H. S. Stern, and D. B. Rubin, 2004: *Bayesian Data Analysis*. 2nd ed. Chapman and Hall, 690 pp.

—, —, —, D. B. Dunson, A. Vehtari, and D. B. Rubin, 2014: *Bayesian Data Analysis*. 3rd ed. Chapman and Hall, 675 pp.

Godart, A., S. Anquetin, and E. Leblois, 2009: Rainfall regimes associated with banded convection in the Cévennes-Vivarais area. *Meteor. Atmos. Phys.*, **103**, 25–34, https://doi.org/10.1007/s00703-008-0326-3.

Liu, J., C. D. Doan, S.-Y. Liong, R. Sanders, A. T. Dao, and T. Fewtrell, 2015: Regional frequency analysis of extreme rainfall events in Jakarta. *Nat. Hazards*, **75**, 1075–1104, https://doi.org/10.1007/s11069-014-1363-5.

Marchi, L., M. Borga, E. Preciso, and E. Gaume, 2010: Characterisation of selected extreme flash floods in Europe and implications for flood risk management. *J. Hydrol.*, **394**, 118–133, https://doi.org/10.1016/j.jhydrol.2010.07.017.

Mélèse, V., J. Blanchet, and G. Molinié, 2018: Uncertainty estimation of Intensity–Duration–Frequency relationships: A regional analysis. *J. Hydrol.*, **558**, 579–591, https://doi.org/10.1016/j.jhydrol.2017.07.054.

—, —, and J.-D. Creutin, 2019: A regional scale-invariant extreme value model of rainfall intensity-duration-area-frequency relationships. *Water Resour. Res.*, **55**, 5539–5558, https://doi.org/10.1029/2018WR024368.
Miniscloux, F., J. D. Creutin, and S. Anquetin, 2001: Geostatistical analysis of orographic rainbands. *J. Appl. Meteor.*, **40**, 1835–1854, https://doi.org/10.1175/1520-0450(2001)040<1835:GAOOR>2.0.CO;2.

Norbiato, D., M. Borga, M. Sangati, and F. Zanon, 2007: Regional frequency analysis of extreme precipitation in the eastern Italian Alps and the August 29, 2003 flash flood. *J. Hydrol.*, **345**, 149–166, https://doi.org/10.1016/j.jhydrol.2007.07.009.

Overeem, A., T. Buishand, I. Holleman, and R. Uijlenhoet, 2010: Extreme value modeling of areal rainfall from weather radar. *Water Resour. Res.*, **46**, W09514, https://doi.org/10.1029/2009WR008517.

Panthou, G., T. Vischel, T. Lebel, G. Quantin, and G. Molinié, 2014: Characterising the space–time structure of rainfall in the Sahel with a view to estimating IDAF curves. *Hydrol. Earth Syst. Sci.*, **18**, 5093–5107, https://doi.org/10.5194/hess-18-5093-2014.

Pauli, F., W. Racugno, and L. Ventura, 2011: Bayesian composite marginal likelihoods. *Stat. Sin.*, **21**, 149–164, http://www.jstor.org/stable/24309266.

Ramos, M. H., J.-D. Creutin, and E. Leblois, 2005: Visualization of storm severity. *J. Hydrol.*, **315**, 295–307, https://doi.org/10.1016/j.jhydrol.2005.04.007.

Ribatet, M., D. Cooley, and A. C. Davison, 2012: Bayesian inference from composite likelihoods, with an application

to spatial extremes. *Stat. Sin.*, **22**, 813–845, http://www.jstor.org/stable/24310036.

Stedinger, J. R., R. M. Vogel, and E. Foufoula-Georgiou, 1993: Frequency analysis of extreme events. *Handbook of Hydrology*, D. R. Maidment, Ed., McGraw-Hill, 18.1–18.66.

Van de Vyver, H., 2015: Bayesian estimation of rainfall intensity–duration–frequency relationships. *J. Hydrol.*, **529**, 1451–1463, https://doi.org/10.1016/j.jhydrol.2015.08.036.

—, and G. R. Demarée, 2010: Construction of Intensity–Duration–Frequency (IDF) curves for precipitation at Lubumbashi, Congo, under the hypothesis of inadequate data. *Hydrol. Sci. J.*, **55**, 555–564, https://doi.org/10.1080/02626661003747390.

Varin, C., 2008: On composite marginal likelihoods. *AStA Adv. Stat. Anal.*, **92**, 1–28, https://doi.org/10.1007/s10182-008-0060-7.

—, and P. Vidoni, 2005: A note on composite likelihood inference and model selection. *Biometrika*, **92**, 519–528, https://doi.org/10.1093/biomet/92.3.519.

Wright, D. B., J. A. Smith, G. Villarini, and M. L. Baeck, 2013: Estimating the frequency of extreme rainfall using weather radar and stochastic storm transposition. *J. Hydrol.*, **488**, 150–165, https://doi.org/10.1016/j.jhydrol.2013.03.003.

Zhou, Z., J. A. Smith, D. B. Wright, M. L. Baeck, L. Yang, and S. Liu, 2019: Storm catalog-based analysis of rainfall heterogeneity and frequency in a complex terrain. *Water Resour. Res.*, **55**, 1871–1889, https://doi.org/10.1029/2018WR023567.