# Modeling Intensity-Duration-Frequency Curves for the Whole Range of Non-Zero Precipitation: A Comparison of Models

Abubakar Haruna, Juliette Blanchet, Anne-catherine Favre

► To cite this version:

Abubakar Haruna, Juliette Blanchet, Anne-catherine Favre. Modeling Intensity-Duration-Frequency Curves for the Whole Range of Non-Zero Precipitation: A Comparison of Models. Water Resources Research, 2023, 59 (6), pp.104044. 10.1029/2022WR033362. hal-04380365

HAL Id: hal-04380365

https://hal.science/hal-04380365v1

Submitted on 18 Jan 2024

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
# Water Resources Research®

## RESEARCH ARTICLE

10.1029/2022WR033362

**Key Points:**

* We build intensity duration frequency (IDF) curves using all the non-zero precipitation data and model the intensities with the extended generalized Pareto distribution

* We consider three approaches to building IDF curves: scale invariance, IDF formulation of Koutsoyiannis et al. (1998), https://doi.org/10.1016/ S0022-1694(98)00097-3 and data-driven method

* We compare and select the best model in a cross-validation framework

**Correspondence to:**

A. Haruna,
abubakar.haruna@univ-grenoble-alpes.fr

**Citation:**

Haruna, A., Blanchet, J., & Favre, A.-C. (2023). Modeling intensity-duration-frequency curves for the whole range of non-zero precipitation: A comparison of models. *Water Resources Research*, 59, e2022WR033362. https://doi.org/10.1029/2022WR033362

Received 31 JUL 2022

Accepted 3 MAY 2023

**Author Contributions:**

**Conceptualization:** Abubakar Haruna, Juliette Blanchet, Anne-Catherine Favre

**Formal analysis:** Abubakar Haruna, Juliette Blanchet, Anne-Catherine Favre

**Methodology:** Abubakar Haruna, Juliette Blanchet, Anne-Catherine Favre

**Software:** Abubakar Haruna

**Supervision:** Juliette Blanchet, Anne-Catherine Favre

**Writing – original draft:** Abubakar Haruna

**Writing – review & editing:** Juliette Blanchet, Anne-Catherine Favre

## Modeling Intensity-Duration-Frequency Curves for the Whole Range of Non-Zero Precipitation: A Comparison of Models

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

Abubakar Haruna¹ ID, Juliette Blanchet² ID, and Anne-Catherine Favre¹

¹University Grenoble Alpes, Grenoble INP, CNRS, IRD, IGE, Grenoble, France, ²University Grenoble Alpes, CNRS, IRD, Grenoble INP, IGE, Grenoble, France

**Abstract** Intensity-duration-frequency (IDF) curves are useful in water resources engineering for the planning and design of hydrological structures. As opposed to the common use of only extreme data to build IDF curves, here, we use all the non-zero rainfall intensities, thereby making efficient use of the available information. We use the extended generalized Pareto distribution to model the distribution of the non-zero rainfall intensities. We consider three commonly used approaches for building IDF curves. The first approach is based on the scale-invariance property of rainfall, the second relies on the general IDF formulation of Koutsoyiannis et al. (1998, https://doi.org/10.1016/S0022-1694(98)00097-3), and the last approach is purely data-driven (Overeem et al., 2008, https://doi.org/10.1016/j.jhydrol.2007.09.044). Using these three approaches, and some extensions around them, we build a total of 10 models for the IDF curves, and then we compare them in a split-sampling cross-validation framework. We consider a total of 81 stations at 10 min resolution in Switzerland. Due to the marked seasonality of rainfall in the study area, we performed a seasonal-based analysis. The results reveal the model based on the data-driven approach as the best model. It is able to correctly model the observed intensities across duration while being reliable and robust. It is also able to reproduce the space and time variability of extreme rainfall across Switzerland.

**Plain Language Summary** Intensity-duration-frequency (IDF) curves are used in water resources engineering for the safe design of hydraulic structures. They summarize the relationship between rainfall intensity, duration, and frequency of occurrence. IDF curves are commonly built using probability distributions that use only extreme data, which can be a drawback in terms of the efficient use of already scarce data. In this study, we use a distribution, the extended generalized Pareto distribution, that uses all the non-zero rainfall intensities to model IDF curves. We build a total of 10 models, which are based on three commonly used approaches and some extensions around them. We compare them in terms of their reliability and robustness. The best IDF model is able to reproduce the known variability of extreme rainfall across Switzerland.

### 1. Introduction

Intensity-duration-frequency (IDF) curves provide the link between precipitation intensity, duration, and non-exceedance frequency (or rather the return period). It is a very common and useful tool in the area of water resources engineering. IDF curves are practically used to infer high return levels of rainfall intensities for the hydrological designs of structures such as sewer lines, culverts, drains, dams, dykes, and so on. They are also used to calibrate/validate stochastic weather generators (Ritschel et al., 2017).

IDF curves are traditionally modeled by first fitting a statistical model, for example, a Gumbel distribution, to extreme data of each duration separately. Second, selected return levels, for example, 2, 5, and 10 years are obtained for each duration from the fitted distribution. And finally, the inferred return levels are linked to duration by some empirical formulation (e.g., Bernard, 1932; Carreteras, 1987; Chow, 1962; Meylan et al., 2012; Sherman, 1931). Although common and easy to implement, there are several drawbacks to this approach. The approach lacks parsimony because several parameters have to be fitted (a set for each return period). It is not robust in the sense that IDF curves are only available for specific return levels. Thus, whenever a new return level is needed, the process has to be repeated. Another major drawback is that uncertainty in the return levels obtained at the initial steps is not taken into account in the last step. Finally, there can be intersections between curves of different return levels that cannot be theoretically justified.

To overcome the outlined limitations of the traditional parametric methods, novel approaches were considered to link the different durations together in IDF curves. In general, in spite of the approach, two choices have to be

© 2023. The Authors.
This is an open access article under
the terms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.
made: the model for the distribution of the rainfall intensities, and the model for the dependence of intensity on duration.

For the choice of the model for the distribution of the rainfall intensities, extreme value distributions are usually considered. For example, in the case of annual maxima series, generalized extreme value (GEV; e.g., Blanchet et al., 2016; Fauer et al., 2021; Innocenti et al., 2017; Jurado et al., 2020; Mélèse et al., 2018; Sane et al., 2018; Ulrich et al., 2020; Van de Vyver, 2018), or a special case of the GEV, that is Gumbel (e.g., Agbazo et al., 2016; Chang et al., 2016; Ghanmi et al., 2016; Yu et al., 2004). For the peaks over threshold (POT), generalized Pareto distribution (GPD) has been used (e.g., Ben-Zvi, 2009; Madsen et al., 1995; Van de Vyver & Demarée, 2010).

Regarding the model for the dependence of intensity on duration, many formulations that are based on different approaches have been proposed in the literature. Here, we identify and focus on three major approaches.

The first approach is based on scale invariance. It has been shown that rainfall exhibits this property within some scales (see Gupta & Waymire, 1990, 1993; Harris et al., 1997; Lima, 1998; Molnar & Burlando, 2008; Over, 1995; Paschalis, 2013; Schertzer & Lovejoy, 1987; Veneziano & Lepore, 2012). This property provides the physical justification for modeling IDF curves, and thus the possibility of inferring return levels of interest across scales. This approach is arguably the most commonly used approach, possibly because of its rich theoretical background, physical basis, and ease of application in regions with scarce availability of sub-daily rainfall series. IDF curves based on this approach can be found in several applications (Blanchet et al., 2016; Burlando & Rosso, 1996; Innocenti et al., 2017; Menabde et al., 1999; Sane et al., 2018; Van de Vyver & Demarée, 2010; Willems, 2000).

The second approach is based on the general formulation of Koutsoyiannis et al. (1998), a generalization of the various empirical formulations for modeling IDF curves. This formulation has the key advantage of being a separable function of return levels and duration. It is also consistent with both probabilistic theories and the physical constraints of scaling across duration. Several applications of this formulation to build IDF curves can be found in the literature (e.g., Blanchet et al., 2016; Fauer et al., 2021; Koutsoyiannis et al., 1998; Roksvåg et al., 2021; Sane et al., 2018; Ulrich et al., 2020; Van de Vyver & Demarée, 2010). This method can be seen as an extension of the scaling approach in which an additional parameter (θ) is added to allow for the curvature of the return levels curves for short durations.

The third approach is based on Overeem et al. (2008) and is called the data-driven approach. Here, the functional relationship (linkage) between IDF model parameters and duration is empirically determined from the data itself. The method involves fitting a parametric model, for example, GEV, to data of each duration. A particular regression model is then fitted for each parameter as a function of duration. As a consequence, the return level of any duration can be inferred from the inverse of the distribution, with parameters obtained from the regression model. This approach imposes neither the assumption/existence of scaling nor the separability condition in the case of the general formulation of Koutsoyiannis et al. (1998). Interestingly, both approaches can be seen as special cases of the data-driven approach with particular functional relationships imposed on the parameters. We note that, although the other two approaches also consider the link between parameter and duration, the specific forms rely on some underlying theoretical hypothesis. However, with the data-driven method, the functional relationship is empirically determined from the data itself, hence its name.

There are also nonparametric approaches, which rather than imposing a parametric model on the intensities, use stochastic rainfall models to estimate the IDF curves (for a brief review, see Langousis & Veneziano, 2007; Veneziano et al., 2007; Tyralis & Langousis, 2019). Here, we focus on the class that uses parametric models for the intensities.

The common use of GEV and GPD as models for the distribution of the rainfall intensities in IDF curves is justified by the practical use of IDF curves, which is to infer high return levels for hydrological designs. However, there is a major drawback in using these distributions. This is mainly due to poor utilization of already scarce data; only one value per block (in GEV distribution), or only values above a given threshold (as in GPD). This is why Langousis et al. (2016) mentioned that compared to the GEV, GPD is more favored for extreme return level estimation since it makes more efficient use of hydrologic information and hence the resulting decrease in estimation uncertainty. Despite that, the delicate issue of threshold remains with the GPD.
As a remedy, Naveau et al. (2016) proposed the extended generalized Pareto distribution (EGPD) to model all the non-zero precipitation intensities. It has the advantage of using all the information present in the sample of non-zero rainfall data, and not only one value per block (like GEV distribution) or only values above a given threshold (as in GPD distribution), thereby reducing uncertainty. It does not require the choice of GPD threshold and has the advantage of being compliant with extreme value theory (EVT) in both the lower and upper tails. Another advantage is its ability to model all non-zero precipitation. This has applications in water resources management for urban water supplies, hydropower, flood forecast, and irrigation systems. As an example, marginal distributions for the entire range of non-zero precipitation are required in stochastic generators (see an example of application in Evin et al., 2018; Viviroli et al., 2022). Simulated precipitation from the generators is used as input to conceptual hydrological models for subsequent flood modeling and risk assessment. Other practical applications are in the evaluation of numerical weather simulations or investigation of the climatology of rainfall events as outlined by Blanchet et al. (2019). On the contrary, the assumption of independence in estimation is more likely to be violated in the case of EGPD, which uses all the information compared to the GEV and GPD. A common approach to reducing the dependence is to apply temporal declustering, for example, by taking one-third of the data as done in Naveau et al. (2016), Le Gall et al. (2022), and Haruna et al. (2022).

The goal of this article is to use the EGPD to build IDF curves for the non-zero precipitation intensities, based on the three outlined approaches, that is, scale-invariance, the general formulation of Koutsoyiannis et al. (1998) and the data-driven approach. Due to the marked seasonality of rainfall in the study area, we consider a seasonally based analysis, with winter (Dec-Jan-Feb), spring (Mar-Apr-May), summer (Jun-Jul-Aug), and autumn (Sep-Oct-Nov). We use a split-sampling cross-validation framework, based on some well-chosen criteria to compare and select the best model. To the best of our knowledge, this study is the first to present an objective comparison of the three outlined IDF curve modeling approaches. Furthermore, while several applications of the EGPD have been presented in the literature (e.g., Blanchet et al., 2019; Evin et al., 2018; Haruna et al., 2022; Le Gall et al., 2022; Rivoire et al., 2021; Tencaliec et al., 2020), our study presents the first attempt to use this distribution in IDF curves development. In addition to the common application of IDF curves for return level estimation, IDF curves based on all the non-zero precipitation data (modeled with EGPD) will allow for a more comprehensive validation of stochastic weather generators. Finally, apart from the traditional IDF curves (see Eicher & Krejci, 1996), as far as we know, this study presents the first attempt to model IDF curves that consider linking parameters with duration over the whole of Switzerland.

The rest of the paper is organized as follows: Section 2 introduces the data and the area under study, and Section 3 presents the EGPD, the models for the IDF curves, and the evaluation framework. Section 4 presents the results and discussion, and finally, Section 5 draws the conclusions and gives relevant perspectives.

## 1.1. Data and Area Under Study

The study area is Switzerland, a small country by size with an area of 41,285 km². Despite its relatively small size, it however presents a complex topography with elevations ranging from 191 to 4,127 m above mean sea level. Around 30% of the area is located above the elevation of 1,500 m above sea level. This results in marked spatial variability both in intensity and occurrence of precipitation.

Point precipitation data from a total of 81 stations, with a minimum record length of 20 years, are available for this study. They are spread across Switzerland and their location is shown in Figure 1. Out of this total, 71 stations belong to the SwissMetNet of the Swiss Federal Office for Meteorology and Climatology (MeteoSwiss), while 10 belong to the canton of Lucerne, a partner network of MeteoSwiss. The precipitation data is measured with a tipping-bucket gauge of 0.1 mm depth resolution at a sampling resolution of 10 min. Most of the stations at high altitudes are shielded from wind and the tipping gauge is heated to account for snow. As a result, heating-related losses can result in up to 24% underestimation compared to the measurements using an electronic weighting system (Savina et al., 2012). The sample data has a variable length ranging from a minimum of 20 years to a maximum of 40 years from 1981 to 2020. The stations are located at elevations ranging from a minimum of 203 m, an average of 952.4 m, and a maximum of 3,294 m.

The climate of Switzerland is influenced by the Alps, the Atlantic Ocean, as well as the Mediterranean sea (Giannakaki & Martius, 2015; Sodemann & Zubler, 2009). As a result, precipitation is characterized by seasonality and spatial variability both in intensity and occurrence (Frei & Schär, 1998; Isotta et al., 2014; Molnar &
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Map of Switzerland showing the location of the 81 stations. The color of the points indicates the length of the precipitation data in years. The background color shows the elevation above sea level in meters.</figcaption></figure>

Burlando, 2008; Sevruk, 1997; Sevruk et al., 1998). Annual precipitation is largest in the Alps and its rims, along the Jura (in the northwest), and in the Ticino (south of the Alps). In these areas, annual sums are generally larger than 2,000 mm. On the other hand, the inner valleys (Rhône and Inn) receive the least annual precipitation (less than 700 mm). Summer is the main season of precipitation all over Switzerland, except in Ticino where autumn is the main season. For all regions, winter receives the least precipitation. In terms of heavy precipitation (see Panziera et al., 2018), the spatial distribution depends on the accumulation duration. For short-duration accumulations (e.g., 1 hr), the heaviest precipitation occurs in summer everywhere. It is largest in Ticino, Jura, and the northern rim, where maximum summer hourly intensities can reach up to 30 mm/hr. For higher accumulations (1 day and greater), Ticino receives the heaviest precipitation, where a maximum of more than 130 mm over 24 hr can be expected in autumn. For the other regions, the heaviest precipitation almost always occurs in summer.

Due to this marked seasonality, we consider a seasonally based analysis. We divide the data into four seasons of 3 months each. Winter includes Dec-Jan-Feb, Spring Mar-Apr-May, Summer Jun-Jul-Aug, while Autumn includes Sep-Oct-Nov. A similar seasonal approach was used in previous studies in this area (see Evin et al., 2018; Fukutome et al., 2015; Haruna et al., 2022; Molnar & Burlando, 2008).

## 2. Methodology

In this section, we start by presenting the distribution for the non-zero precipitation intensities, we then present the various IDF models, and finally, the inference strategy to estimate the parameters. A flowchart for the framework is shown in Figure 2.

### 2.1. Marginal Distribution for Non-Zero Precipitation Intensities

We define the random variable *I* to represent non-zero rainfall intensities. We intend to find a distribution that will allow us to model the entire distribution of *I*, that is, both its small, medium, and extreme values. A natural choice would be the Gamma distribution, but it has been reported to underestimate extreme precipitation (Katz et al., 2002).

When looking at only the extremes, the EVT (Coles, 2001) tells us that, the probability that the excesses of *I* above a high threshold *u* are less than or equal to a large value *i* can be approximated by the GPD given as:

$$ \mathbb{P}(I - u \le i | I > u) \rightarrow H_{\xi}\left(\frac{i}{\sigma}\right) \quad (1) $$

with the cumulative distribution function (CDF)
<figure><img src="image_4.png" /><figcaption>Figure 2. Flowchart of the framework for the study.</figcaption></figure>

$$H_{\xi}\left(\frac{i}{\sigma}\right)=\begin{cases}1-\left(1+\xi \frac{i}{\sigma}\right)^{-1 / \xi} & \text { if } \xi \neq 0 \\ 1-\exp \left(-\frac{i}{\sigma}\right) & \text { if } \xi=0\end{cases} \quad (2)$$

where $a_+ = \max(a, 0)$, $\sigma > 0$ is the scale parameter, and $-\infty < \xi < +\infty$ is the GP shape parameter that controls the upper tail of the distribution. The bounded case (short-tailed) is when $\xi < 0$, the exponential case (light-tailed) when $\xi = 0$, and the unbounded case (heavy-tailed) when $\xi > 0$.

The GPD has a key advantage over the GEV since it makes use of all the observations above the threshold $u$, compared to the GEV case where only the largest observation within a block, usually a season or year, is retained. However, some major issues remain. First, there is yet no unified method for the choice of the optimum threshold $u$ over which excesses are GP distributed. Second, the observations that are below the chosen threshold $u$, although precious and not easy to come by, are discarded and entirely not utilized. Finally, the question of modeling the whole range of non-zero precipitation which has applications in water resource management remains.

Many parametric, semi-parametric, and non-parametric models have been proposed in the literature as extensions to the GPD to model the full distribution of $I$, and not only its upper tail as given in Equation 1 (see review in
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Density of the extended generalized Pareto distribution (EGPD) for $G(v) = v^k$ for $\sigma = 1$, $\xi = 0.5$, and $\kappa = 0.8, 1$, and 2 (adapted from Naveau et al. (2016)). The density of gamma distribution with both shape and scale = 1.4 is shown in blue color. The right panel zooms on the right tail. Indeed the gamma tail decays much faster while having a similar shape in the bulk of the distribution for the case of $\kappa = 2$. The case where $\kappa = 1$ is the exact GP case.</figcaption></figure>

Scarrot & MacDonald, 2012). In particular, Naveau et al. (2016) proposed a family of models they called the EGPDs that are defined as:

$$ \mathbb{P}(I \le i) = G\left[H_{\xi}\left(\frac{i}{\sigma}\right)\right], \qquad (3) $$

where G is a continuous CDF that is also defined on [0,1], with constraints given in Naveau et al. (2016) to ensure that the distribution of I remains GP in its upper tail, a gamma-like lower tail, and a smooth transition in between. Four parametric families of G were proposed by the authors, the simplest of which is defined as $G(v) = v^k$. The CDF of this model is defined as:

$$ \mathbb{P}(I \le i) = \left[ H_{\xi}\left(\frac{i}{\sigma}\right) \right]^k \qquad (4) $$

The model thus has three parameters. $k > 0$ controls the lower tail, $\sigma > 0$ is the scale parameter, and $\xi \ge 0$ controls the upper tail. $H_{\xi}\left(\frac{i}{\sigma}\right)$ is as defined in Equation 2. We exclude the bounded case ($\xi < 0$) since we are dealing with precipitation. We thus assume precipitation to be either light-tailed ($\xi = 0$) or heavy-tailed ($\xi > 0$). The density of the model, alongside that of the Gamma distribution, is illustrated in Figure 3. With just one additional parameter, $\kappa$, the distribution is parsimonious and able to adequately model the full range of non-zero precipitation (see applications in Evin et al., 2018; Haruna et al., 2022; Le Gall et al., 2022). We thus use this model in the rest of the paper.

## 2.2. IDF Models

We define the random variable $I_d$ as the average non-zero precipitation intensity over the duration d. It is described by the CDF, $F_d(i)$, such that $F_d(i) = \mathbb{P}(I_d < i)$. The exceedance frequency is defined as $p_d(i) = 1 - F_d(i)$. The return period of any non-zero intensity i, as a function of $p_d$ is given by $T(I_d \ge i) = \frac{1}{p_d \times \delta_d}$, with $\delta_d$ the average number of non-zero precipitation intensities of duration d per year. We estimate $\delta_d$ based on the long-term average of the non-zero precipitation intensities per year. Consequently, the T-year return level over duration d, $i(T, d)$, is defined as the $(1 - \frac{1}{T \times \delta_d})$ quantile of $F_d$.

Accordingly, IDF relationship is a mathematical function $(T, d) \mapsto i(T, d)$ that relates non-zero rainfall intensity i with its duration d, and the frequency of exceedance $p_d$ (or rather the return level T). In this article, the CDF of $I_d$, $F_d(i)$ is defined by the EGPD presented in Section 3.1. All the different formulations considered here simply differ by how they define this mathematical relationship between i, T, and d, while taking $F_d(i)$ as an EGPD model.
In the following subsections, we present the different IDF-EGPD models based on the three outlined approaches, that is, scale-invariance, the general formulation of Koutsoyiannis et al. (1998), and data-driven approaches. For sake of simplicity, we drop the “EGPD” term and simply refer to the IDF models as IDF<sub>modelname</sub>, where the subscript “modelname” refers to the approach used to build the model.

For all the models, the IDF curve, corresponding to the $(1 - \frac{1}{T \times \delta_d})$ quantile of the EGPD is defined in Equation 5 as:

$$i(T, d) = \frac{\sigma_d}{\xi_d} \left\{ \left[ 1 - \left[ 1 - \frac{1}{T \times \delta_d} \right]^{\frac{1}{\kappa_d}} \right]^{-\xi_d} - 1 \right\} \quad (5)$$

Where $\kappa_d > 0$, $\sigma_d > 0$, and $\xi_d \ge 0$ are the three EGPD parameters for the duration $d$, $T$ is the return period in years, and $\delta_d$ is the average number of non-zero precipitation intensities per year for the duration $d$. The choice of the model determines whether each of the three parameters; $\kappa$, $\sigma$, and $\xi$ varies with $d$ or not, and the form of the relationship.

We consider 13 durations, that is, $d = 30$ min, 40 min, 1, 2, 3, 6, 10, 12, 16, 18, 24, 48, and 72 hr. We consider durations up to 72 hr (3 days) because according to Froidevaux et al. (2015), precipitation accumulations from 0 to 3 days before an event are the most relevant for triggering floods of high magnitudes in Switzerland. The intermediate durations are meant to ensure a good spread in a logarithmic scale and for a later comparison of our quantiles with those provided by MeteoSwiss (not in this paper). We use a fixed window to aggregate the data from the gauge resolution of 10 min, to the various durations. For instance, the 24 hr intensities correspond to amounts accumulated from 00h00 to 24h00 UTC of every day, divided by 24. We chose this time window because, in the study area, convective events are mostly in the evening. This will allow us to differentiate between the first and last 12 hr starting from 00 hr. The choice of a fixed window over a moving window will likely result in omitting the highest intensities in each duration. However, since we are using all the non-zero precipitation intensities, using a moving window will result in significant dependence in the time series.

### 2.2.1. Scaling IDF Model

Scale invariance in the strict sense of Gupta and Waymire (1990) refers to the property where the probability distribution of $I_d$ can be inferred from the distribution of $I_{d_0}$ at the reference duration $d_0$ through:

$$I_d \stackrel{dist}{=} C_\lambda I_{d_0}, \quad (6)$$

where the parameter $C_\lambda$ determines the type of scaling; simple-scaling or multi-scaling. For our case, the reference duration $d_0$ is taken as 1 hr.

A weaker assumption, the so-called “wide sense scaling” (Gupta & Waymire, 1990), is when the scaling is in the moments according to:

$$\mathbb{E}[I_d^q] = \left(\frac{d}{d_0}\right)^{-k(q)} \mathbb{E}[I_{d_0}^q], \quad (7)$$

where $q$ is the order of the moment, $k(q)$ is called the moment scaling function, and $d_0$ is the reference duration. Moment scaling analysis as described by Gupta and Waymire (1990) is used to determine the type of scaling.

When the parameter in Equation 6, $C_\lambda = (\frac{d}{d_0})^{-H}$, that is, a scalar that depends only on the ratio of the scales, we have “strict sense simple-scaling.” This is expressed in Equation 8.

$$I_d \stackrel{dist}{=} \left(\frac{d}{d_0}\right)^{-H} I_{d_0} \quad (8)$$

where $0 < H < 1$ is the scaling exponent. The other variables retain their meanings.

Wide sense simple-scaling is when the moment scaling function in Equation 7 is linear in $q$, that is, $k(q) = Hq$, as expressed in Equation 9.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Illustration of (a) Break in scaling of the $\sigma$ parameter in winter at a station, Robbia in Graubünden. The points colored black are the estimated $\sigma$ for each duration separately. The models in blue and red are given in Equations 10 and 11, respectively. (b) Dependence of $\xi$ on duration in summer at a station in Zurich. The black colored points are the estimated $\xi$ for each duration separately. The lines are the fitted linear models whose equations are given in Appendix B. The broken lines represent the 95% confidence interval for the fitted models. The plots of the complete extended generalized Pareto distribution parameters with duration at these two stations are also given in Appendix B2 (see Figures B1–B2B3B4B5B6) for the four seasons.</figcaption></figure>

$$ \mathbb{E}[I_d^q] = \left(\frac{d}{d_0}\right)^{-Hq} \mathbb{E}[I_{d_0}^q] \qquad (9) $$

It can be shown that, under the strict sense simple-scaling, only one parameter of the EGPD is scaling, which is $\sigma$, whereas $\kappa$ and $\xi$ are independent of duration. For the rest of the paper, we drop the term “strict sense,” and simply use “simple scaling” for convenience. Accordingly, the simple-scaling EGPD model, IDF$_{ss}$, is defined such that: $\kappa_d = \kappa_{d_0}$, $\xi_d = \xi_{d_0}$, and $\sigma_d$ is a power law given as:

$$ \sigma_d = \left(\frac{d}{d_0}\right)^{-H} \sigma_{d_0} \qquad (10) $$

where $\kappa_{d_0} > 0$, $\sigma_{d_0} > 0$, and $\xi_{d_0} \ge 0$, are the parameters of the reference duration $d_0 = 1$ hr, and $0 < H < 1$ is the scaling exponent. The inference method for the parameters is described in Section 3.3.

An important issue is the existence of multiple scaling regimes in precipitation. This means that different scaling behaviors (scaling exponents) exist for different ranges of duration. IDF curves have to be modeled considering the existence of this change in scaling (e.g., Bougadis & Adamowski, 2006; Courty et al., 2019; Yu et al., 2004). To illustrate this behavior, we consider the case of a station at Robbia in Graubünden in winter. We fit the EGPD to the data of each of the 13 durations separately by maximum likelihood. We then inspected how the estimated scale parameter ($\sigma$) varies with duration. The result is given in Figure 4a. Here, a single power law (log-log given in Equation 10) in blue is not enough to explain the scaling. To account for this break in the scaling relationship, we define the two-regime simple-scaling EGPD IDF model, IDF$_{ss\_TR}$ as:

$$ \sigma_d = \begin{cases} \left(\frac{d}{d_0}\right)^{-H_1} \sigma_{d_0} & \text{if } d \le K, \\ \left(\frac{d}{d_0}\right)^{-H_2} \sigma_{d_0} \times K^{H_2-H_1} & \text{if } d > K \end{cases} \qquad (11) $$

An illustration of this equation is in Figure 4a (log-log TR in red). Where K is the duration of the scaling break, $\sigma_{d_0} > 0$ is the scale parameter for the reference duration $d_0 = 1$ and $\sigma_d$ is continuous in $d = K$. The other parameters, $0 < H_1 < 1$ and $0 < H_2 < 1$ are the scaling exponents of the first and second regimes. The other two EGPD
parameters $κ_d = κ_{d0}$, $ξ_d = ξ_{d0}$, remain independent of durations. Hence for this model, a total of six parameters have to be estimated, that is, $κ_{d0}$, $σ_{d0}$, $ξ_{d0}$, $H_1$, $H_2$ and K.

Finally, although the simple-scaling EGPD model imposes a constraint on the linkage of $ξ$ with $d$, that is, $ξ_d = ξ_{d0}$, we however notice some of the stations to show an apparent variation of $ξ$ with respect to $d$. For illustration, we consider a station in Zurich in summer. We fit EGPD to the data of each of the 13 durations separately by maximum likelihood. We then inspected how the estimated shape parameter ($ξ$) varied with duration and modeled the relationship through a linear-log form as expressed in Equation 12. Figure 4b illustrates this and how the linear-log model fits the points correctly:

$$ξ_d = a_ξ + b_ξ \log(d), \qquad (12)$$

where $a_ξ$ and $b_ξ$ are the intercepts and slopes, respectively. This leads to two additional IDF models, with $ξ = f(d)$, namely:

* IDF$_{ss_ξ(d)}$: an extension of the basic simple-scaling model IDF$_{ss}$, to allow $ξ$ to vary with $d$ according to Equation 12.

* IDF$_{ss_TR_ξ(d)}$: an extension of the two-regime simple-scaling model IDF$_{ss_TR}$, to allow $ξ$ to vary with $d$ according to Equation 12.

### 2.2.2. General IDF Formulation

Koutsoyiannis et al. (1998) proposed a general formulation for the different traditional formulations of the IDF curves in the literature. He showed that all of them can be simplified into the form:

$$i(T, d) = \frac{a(T)}{b(d)} \qquad (13)$$

where $b(d) = (d + θ)^H$. The parameter $θ > 0$ is the duration offset, and $0 < H < 1$ is the duration exponent. $a(T)$ is the $(1 - 1/T)$ quantile of the re-scaled intensities $I_d b(d)$. $a(T)$ is independent of $d$ and completely determined by the statistical model considered for $I_d$, in our case, the EGPD.

This formulation has the key advantage of being a separable function of return levels $a(T)$, and duration $b(d)$ that is consistent with both probabilistic theories and the physical constraints of scaling across duration. Menabde et al. (1999) showed that this formulation is the same as the scale-invariant model if $θ$ is set to 0.

When applied to the EGPD, IDF$_{kouts}$ is defined such that: $κ_d = κ_{d0}$, $σ_d = (d + θ)^{-H_1} σ_{d0}$, $ξ_d = ξ_{d0}$. Five parameters, $κ_{d0} > 0$, $σ_{d0} > 0$, $ξ_{d0} ≥ 0$, $θ > 0$, and $0 < H < 1$ have to be inferred. The reference duration here is taken as $d_0 = 1 - θ$.

Following the same arguments discussed in Section 3.2.1 regarding the existence of a break in the scaling relationship, and the dependence of $ξ$ with $d$, we propose three extensions to this model:

* IDF$_{kouts_TR}$: Allowing for a break in the scaling regime. This model is defined as:

$$\sigma_d = \begin{cases} \left(\frac{d + \theta}{d_0 + \theta}\right)^{-H_1} \sigma_{d_0} & \text{if } d \le K \\ \left(\frac{d}{d_0}\right)^{-H_2} \sigma_{d_0} \times K^{H_2 - H_1} & \text{if } d > K \end{cases} \qquad (14)$$

where K is the duration of the scaling break, $σ_{d0} > 0$ is the scale parameter for the reference duration $d_0 = 1 - θ$ and $σ_d$ is continuous in $d = K$. The other parameters, $0 < H_1 < 1$ and $0 < H_2 < 1$ are the scaling exponents of the first and second regimes. The other two EGPD parameters $κ_d = κ_{d0}$, $ξ_d = ξ_{d0}$, remain independent of durations. Hence for this model, a total of six parameters have to be estimated, that is, $κ_{d0}$, $σ_{d0}$, $ξ_{d0}$, $H_1$, $H_2$, and K.

* IDF$_{kouts_ξ(d)}$: an extension of the basic model IDF$_{kouts}$, to allow $ξ$ to vary with $d$ according to Equation 12.

* IDF$_{kouts_TR_ξ(d)}$: an extension of the two-regime model IDF$_{kouts_TR}$ to allow $ξ$ to vary with $d$ according to Equation 12.
<table><caption>**Table 1**
Summary of the IDF Models That Are Compared in This Study</caption><thead><tr><th></th><th>Model</th><th>No. of parameters</th><th>Name of approach</th><th>Inference method</th><th>Ref. section</th></tr></thead><tbody><tr><td>1</td><td>IDF<sub>ss</sub></td><td>4</td><td>Simple-scaling</td><td>Global MLE<sup>a</sup></td><td>3.2.1</td></tr><tr><td>2</td><td>IDF<sub>ss_TR</sub></td><td>6</td><td>Simple-scaling</td><td>" "</td><td>3.2.1</td></tr><tr><td>3</td><td>IDF<sub>ss_ξ(d)</sub></td><td>5</td><td>Extension of Simple-scaling</td><td>" "</td><td>3.2.1</td></tr><tr><td>4</td><td>IDF<sub>ss_TR_ξ(d)</sub></td><td>7</td><td>Extension of Simple-scaling</td><td>" "</td><td>3.2.1</td></tr><tr><td>5</td><td>IDF<sub>kouts</sub></td><td>5</td><td>Koutsoyiannis et al. (1998)</td><td>" "</td><td>3.2.2</td></tr><tr><td>6</td><td>IDF<sub>kouts_TR</sub></td><td>7</td><td>Ext. of Koutsoyiannis et al. (1998)</td><td>" "</td><td>3.2.2</td></tr><tr><td>7</td><td>IDF<sub>kouts_ξ(d)</sub></td><td>6</td><td>Ext. of Koutsoyiannis et al. (1998)</td><td>" "</td><td>3.2.2</td></tr><tr><td>8</td><td>IDF<sub>kouts_TR_ξ(d)</sub></td><td>8</td><td>Ext. of Koutsoyiannis et al. (1998)</td><td>" "</td><td>3.2.2</td></tr><tr><td>9</td><td>IDF<sub>DD<sub>global</sub></sub></td><td>10</td><td>Data-driven</td><td>" "</td><td>3.2.3</td></tr><tr><td>10</td><td>IDF<sub>DD<sub>local</sub></sub></td><td>10</td><td>Data-driven</td><td>Two-step(MLE + regression)</td><td>3.2.3</td></tr></tbody></table>

### 2.2.3. Data-Driven IDF Model

The scaling theory and the specific form of Equation 13 impose particular functions for the relation between the scale parameter, $\sigma$ of the EGPD with respect to duration, $d$. However, in the case of the data-driven models, the expression of the relationship for each of the three EGPD parameters is empirically determined by the data itself. To guide our choice of the appropriate functional relationship, we inspected how each locally estimated EGPD parameter varies with duration. Figure 4 gives an example for the $\sigma$ and $\xi$ parameters at two stations. We finally settled on the following functions to model the three parameters with respect to duration:

$$ \kappa_d = \begin{cases} \exp[a_\kappa + b_{1,\kappa}\log(d)] & \text{if } d \le K_\kappa \\ \exp[a_\kappa + b_{2,\kappa}\log(d) + (b_{1,\kappa} - b_{2,\kappa})\log(K_\kappa)] & \text{if } d > K_\kappa \end{cases} \quad (15) $$

$$ \sigma_d = \begin{cases} \exp[a_\sigma + b_{1,\sigma}\log(d)] & \text{if } d \le K_\sigma \\ \exp[a_\sigma + b_{2,\sigma}\log(d) + (b_{1,\sigma} - b_{2,\sigma})\log(K_\sigma)] & \text{if } d > K_\sigma \end{cases} \quad (16) $$

For the first two parameters, $\kappa$ and $\sigma$, the function is a continuous two-linear piece-wise model in log space. $K_*$ is the duration of the break-point ($\sigma$ continuous for $d = K_*$). $a_*$, $b_{1,\*}$, and $b_{2,\*}$ are the intercepts and slopes of the first and second lines, respectively. In the case of $\xi$, the function was already given in Equation 12.

Note that, by keeping $\kappa$ and $\sigma$ independent of duration, and using either $\sigma_d = (\frac{d}{d_0})^{-H} \sigma_{d_0}$ or $\sigma_d = (\frac{d+\theta}{d_0+\theta})^{-H} \sigma_{d_0}$, the simple-scaling or the general formulations of Koutsoyiannis et al. (1998) presented in Sections 3.2.1 and 3.2.2 respectively can be obtained from this data-driven approach.

We consider two IDF models in this class, $IDF_{DD_{global}}$ and $IDF_{DD_{local}}$, both impose the same type of functional relationships (Equations 12, 15, and 16), but simply differ in the way the regression parameters are estimated. The inference strategy is explained in detail in Section 3.3.

The different models compared in this study are summarized in Table 1. For the list of some key variables and their meanings, refer to Table A1.

### 2.3. Inference

In this section, we describe the inference methods to estimate the parameters of the ten (10) IDF models presented in Table 1.

<footnote>aMaximum Likelihood Estimation.</footnote>
As a prerequisite to an objective comparison, the same estimation strategy has to be employed for all the models. For this reason, we use a global maximum likelihood estimation for all the models (as done by Blanchet et al. (2016)), which we describe in the next paragraph. An exception is the case of only one model, the IDF<sub>DD<sub>local</sub></sub>, which involves a two-step method. In this case, first, we fit, for a given station, the EGPD on the data of each duration separately (using maximum likelihood estimation). Second, we fit for each fitted parameter, the chosen regression model as a function of duration. The parameters in Equations 15 and 16 are estimated by segmented regression (i.e., regression model with break-points), while those of Equation 12 by ordinary least squares. The segmented regression we use here is based on the algorithm of Muggeo (2003), and with the functions that are implemented in the segmented package. Details of the algorithm can be found in Muggeo (2003).

We now come back to the global maximum likelihood for the other nine (9) IDF models. This method involves pooling, for each station, all the data from the 13 durations to estimate the model parameters. The duration *d* is used as a covariate. We note here that by pooling all the data, we made the hypothesis of independence between the intensities of the different time steps and different durations (since we use independence likelihood as in Blanchet et al. (2016)). This hypothesis is difficult to justify given that we consider all the non-zero intensities. Note that by taking fixed windows rather than moving windows, we have reduced the dependence between the different time steps. Previous studies considered similar approaches to reduce temporal dependence between different time steps, by taking a fraction of the data (Haruna et al., 2022; Le Gall et al., 2022; Naveau et al., 2016)). The dependence between the intensities of different duration, for example, intensities of 1 hr versus those of 2 hr, however, remains. While it might be possible to account for this dependence using different strategies proposed in the literature (e.g., Davison et al., 2012; Sebille et al., 2017), we believe this requires a complete task of its own and will likely complicate the optimization of the models. This is in addition to the possibility of performance deterioration due to dependence structure misspecification. We take this as a limitation of the present study and a topic for future research.

The log-likelihood (*ll*) that is maximized here (given in Equation 17) takes left censoring into account. The importance of using left censoring in fitting rainfall data by maximum likelihood has been pointed out by Naveau et al. (2016), and they showed that better performance is obtained by taking it into account.

$$ll_{EGPD}(\kappa_d, \sigma_d, \xi_d) = ll_{censored}(\kappa_d, \sigma_d, \xi_d) + ll_{uncensored}(\kappa_d, \sigma_d, \xi_d), \quad (17)$$

where *ll*<sub>censored</sub> and *ll*<sub>uncensored</sub> are the contributions of the censored and uncensored data, given in Equations 18 and 19, respectively, as

$$ll_{censored}(\kappa_d, \sigma_d, \xi_d) = \sum_d \sum_{j:i_d \le c_d} \kappa_d \log \left[ 1 - \left( 1 + \frac{\xi_d c_d}{\sigma_d} \right)^{-\frac{1}{\xi_d}} \right], \quad (18)$$

$$ll_{uncensored}(\kappa_d, \sigma_d, \xi_d) = \sum_d \sum_{j:i_d \ge c_d} \log \kappa_d - \sum_d \sum_{j:i_d \ge c_d} \log \sigma_d - \sum_d \sum_{j:i_d \ge c_d} \left[ 1 + \frac{\xi_d i_{d,j}}{\sigma_d} \right]^{\left[ 1 + \frac{1}{\xi_d} \right]} + \\ \sum_d \sum_{j:i_d \ge c_d} \left[ 1 - \left( 1 + \frac{\xi_d i_{d,j}}{\sigma_d} \right)^{-\frac{1}{\xi_d}} \right]^{[\kappa_d - 1]}, \quad (19)$$

where *d* ranges over the 13 durations and *j* ranges over the number of time steps in the data of duration *d*. *κ*<sub>*d*</sub> > 0, *σ*<sub>*d*</sub> > 0, *ξ*<sub>*d*</sub> ≥ 0 are the EGPD parameters for duration *d*. *i*<sub>*d,j*</sub> is the precipitation intensity for the duration *d* and time step *j*. *c*<sub>*d*</sub> ≥ 0 is the left censoring threshold applied to the data of duration *d*. Many authors have taken this censoring approach into account but they usually take a uniform threshold value for all the stations (e.g., Tencaliec et al. (2020) used 2 mm for daily rainfall). Here, we did not find the use of a common threshold over the 81 stations sufficient. We had to select, for each station and duration, the lower threshold *c* that minimizes the Normalized Root Mean Square Error (NRMSE) of Equation 20 in Section 3.5.

In both Equations 18 and 19, the choice of the IDF model specifies the function linking the EGPD parameters to duration. For instance, in the case of the simple-scaling, IDF<sub>DD<sub>ss</sub></sub>, κ<sub>d</sub> = κ<sub>d0</sub> σ<sub>d</sub> = (d/d<sub>0</sub>)<sup>-H</sup> σ<sub>d0</sub> ξ<sub>d</sub> = ξ<sub>d0</sub>.
For the cases where two linear models are fitted (see Equations 11, 12, 14–16), we use segmented regression to estimate the regression parameters. We then use the fitted parameters as initial values in the optimization of the likelihood function.

We note here that other estimation methods were used in other studies for the simple-scaling approach and the general IDF formulation of Koutsoyiannis et al. (1998), besides the global MLE. In the case of the simple-scaling models, a two-step procedure, where the scaling exponent in Equation 8 is first obtained through moment scaling analysis, then all re-scaled intensities from all the durations are pooled to fit the IDF model using MLE (see Innocenti et al., 2017; Nhat et al., 2008; Panthou et al., 2014). For the general IDF formulation of Koutsoyiannis et al. (1998), the authors proposed two different estimation strategies; the so-called “robust estimation” and the “one-step least square method”. The robust estimation is a two-step procedure that involves the estimation of the parameters of $b(d)$, and then those of $a(T)$ (see Equation 13), through the minimization of the Kruskal-Wallis statistic. The one-step least square method involves the joint estimation of all the parameters of Equation 13 that minimizes the squared error of the observed and modeled quantiles from the IDF model. Despite this, we only use the global MLE method to objectively compare the models based on the same inference strategy. In addition, it has the advantage of being a one-step estimation procedure and is better suited than least squares for non-Gaussian distribution estimation.

## 2.4. Uncertainty Estimation

We employ a block bootstrap approach for the uncertainty estimation in all the models considered in this study. The principle of the block bootstrap involves sampling with replacement, all the data contained in a block of a given size *B*, *R* number of times. Here, for computational reasons, we use *R* = 500. We then use the percentile method to estimate the 95% confidence interval (CI). The choice of the appropriate block size *B* is a delicate issue in the scientific community, but one common way is to choose *B* large enough to ensure that the temporal dependence in the data is maintained. A block size *B* too small will underestimate the uncertainty. Here, we checked the autocorrelation in the seasonal data (result not shown), and saw that it does not decrease after a lag of more than 1 week. Hence we choose *B* = 2 weeks. To maintain the dependence between data of different durations *d*, we always ensured that all the data of the different durations *d* contained in the same block *B* = 2 weeks are sampled together. A summary of the block bootstrapping is summarized in the steps below:

For each station *s* and season:

1. Aggregate the data to intensities of the 13 duration *d*, store the data in a matrix of *n* by *d*. Where *n* is the number of observations and *d* = 13. We call this matrix, $M_{orig}$.

2. Sample with replacement, $\frac{n}{B}$ blocks from the row of matrix $M_{orig}$ in the previous step to form the bootstrap matrix $M_{boot}$. Both $M_{orig}$ and $M_{boot}$ have the same dimensions. By sampling from the rows of $M_{orig}$, we keep the data of the different durations *d* together, and hence the dependence structure.

3. Fit the IDF model on the data in $M_{boot}$ and estimate the intended return level.

4. Repeat steps 2 to 3 a total of 500 times ($R = 500$) to obtain the bootstrap distribution of the return level.

5. Obtain the 95% CI of the intended return level by percentile method. This is done by taking the empirical 0.025 and 0.975 quantiles of the bootstrap distribution of the return level obtained in step 4.

We note here that to ensure all the models use the same bootstrap matrix $M_{boot}$, we use the “set.seed” in to keep track of the random number generation in the sampling of step 2. Examples of the CI obtained with this approach are later given in Figure 10.

## 2.5. Evaluation Framework

We evaluate the performance of the models in two aspects. First, in calibration, that is how well a given model predicts the data that was used in training it. Second, we evaluate their predictive performance in a cross-validation framework.
## 2.5.1. Calibration

We use two criteria to evaluate the performance of the models in calibration (in-sample performance). They are itemized below:

* NRMSE: The normalization, which here is done by the mean, allows the comparison of intensities of different duration across different stations. For each station s, and duration d, we compute the NRMSE over the non-zero precipitation intensities as:

$$ NRMSE_s(d) = \frac{\left\{ \frac{1}{n_s(d)} \sum_{j=1}^{n_s(d)} (r_{s,T_j}(d) - \hat{r}_{s,T_j}(d))^2 \right\}^{1/2}}{r_s(d)} \quad (20) $$

where NRMSE_s(d) denotes the score computed at station s, and duration d, n_s(d) is the number of non-zero precipitation intensities for duration d, r_{s,T_j}(d) is the empirical quantile with return period T_j = (n_s(d) + 1)/(j × δ_d), δ_d the average number of non-zero precipitations of duration d per year, r_{s,T_j}(d) is the corresponding T_j-year return level estimated from the fitted model. The denominator is the average rainfall at site s and duration d, given as (1/n_s(d)) ∑_{j=1}^{n_s(d)} r_{s,T_j}(d). The best model according to this criteria is the one with the lowest NRMSE_s(d).

* Akaike Information Criteria (AIC) (Akaike, 1974): This information criterion rewards the goodness of fit, as measured by the likelihood, but also penalizes the additional number of parameters to be estimated. It is computed as:

$$ AIC = -2 \log(L) + 2p \quad (21) $$

where L is the maximized likelihood from Equation 17 and p is the number of parameters to be estimated (see Table 1). The lower the AIC, the better the model.

## 2.5.2. Cross-Validation

We follow a split sampling procedure in a cross-validation framework. For each station s, we divide the 10 min precipitation intensities into two equal sub-samples of the same length but on different years that are randomly chosen. We then aggregate the data into intensities of various duration, d = 30 min, 40 min, 1,2, 3, 6, 10, 12, 16, 18, 24, 48, and 72 hr. Then we fit each of the 10 IDF models.

We then evaluate the performance of the models fitted on sub-sample 1 on the observations in sub-sample 2 and vice versa, by computing the relevant cross-validation criteria. We repeat this procedure 10 times.

In the following, we present three different criteria we use to evaluate the models. These criteria have been used in several studies to evaluate and compare competing models (see Blanchet et al., 2015; Evin et al., 2016; Garavaglia et al., 2011; Haruna et al., 2022; Renard et al., 2013).

* The Robustness criteria, SPAN, measures the stability of the estimate of a high return level when the training data is changed. It is computed as:

$$ SPAN_{s,T}(d) = \frac{2|\hat{r}_{s,T}^{(1)}(d) - \hat{r}_{s,T}^{(2)}(d)|}{(\hat{r}_{s,T}^{(1)}(d) + \hat{r}_{s,T}^{(2)}(d))} \quad (22) $$

where $\hat{r}_{s,T}^{(1)}(d)$ and $\hat{r}_{s,T}^{(2)}(d)$ are the T-year return levels estimated from sub-samples 1 and 2, respectively, at station s and duration d. A SPAN of 0.5 means that the absolute difference between the two return levels is half of their average.
For each duration *d*, a regional score over all the *N* stations (*N* = 81) is computed as $SPAN_{reg,T}(d) = 1 - \frac{1}{N} \sum_{s=1}^{N} SPAN_{s,T}(d)$ and a perfect model in terms of robustness according to this criteria should have $SPAN_{reg,T}(d) = 1$.

* The reliability of the model fitted on sub-sample 1 in predicting the maxima in sub-sample 2 and vice versa is measured by the FF criteria:

$$ FF_s^{(12)}(d) = [\hat{F}_s^{(1)}(d)(\max_s^{(2)}(d))]^{n_s^{(2)}(d)} \quad (23) $$

where $FF_s^{(12)}(d)$ is the cross-validation criteria computed at station *s*, and duration *d*, by predicting the probability of the maximum value in sub-sample 2, of sample size $n_s^{(2)}(d)$ using the model $\hat{F}_s^{(1)}(d)$ fitted on the sub-sample 1. $FF_s^{(21)}(d)$ is computed symmetrically.

For a given duration, Renard et al. (2013) and Blanchet et al. (2015) showed that each $FF_s^{(12)}(d)$ at a station should be a realization of a uniform distribution. So the difference in the area, *diff* between a theoretical uniform distribution and that of the *N* set of $FF_s^{(12)}(d)$ (computed over the *N* stations), should be close to 0. $FF_{reg}(d)$ at the regional scale, given as $1 - diff$, should therefore take a value of 1 for a reliable model and 0 for a completely unreliable model; the lower the value, the less reliable the model is.

* The reliability/accuracy of the model in predicting the entire observations in cross-validation is measured by the NRMSE_CV.

$$ NRMSE_CV_s^{(12)}(d) = \frac{\left\{ \frac{1}{n_s^{(2)}(d)} \sum_{j=1}^{n_s^{(2)}(d)} \left( r_{s,T_j}^{(2)}(d) - \hat{r}_{s,T_j}^{(1)}(d) \right)^2 \right\}^{1/2}}{\hat{r}_s^{(2)}(d)} \quad (24) $$

where $NRMSE_CV_s^{(12)}(d)$ is the score computed at station *s*, and duration *d*, $n_s^{(2)}(d)$ is the sample size, $r_{s,T_j}^{(2)}(d)$ is the empirical quantile with return period $T_j$ in sub-sample 2, $\hat{r}_{s,T_j}^{(1)}(d)$ is the corresponding $T_j$ return level estimated from $\hat{F}_s^{(1)}(d)$. The denominator is the average daily rainfall in sub-sample 2 at site *s* given as $\frac{1}{n_s^{(2)}(d)} \sum_{j=1}^{n_s^{(2)}(d)} r_{s,T_j}^{(2)}(d)$.

Finally, for each duration *d*, the regional score computed over the *N* stations is given as: $NRMSE_CV_{reg}^{(12)}(d) = 1 - \frac{1}{N} \sum_{s=1}^{N} NRMSE_CV_s^{(12)}(d)$. $NRMSE_CV_{reg}^{(21)}(d)$ is computed in similar way. $NRMSE_CV_{reg} = 1$ means a perfect model and indicates that there is a complete agreement between all the empirical and theoretical quantiles for all return periods. The closer the value is to 1, the more accurate the model is.

# 3. Results and Discussion

We present the results in the following order: first, we investigate the appropriateness of the EGPD to fit the data of each duration. Then we present the results of the comparison of the IDF models in calibration, and then in cross-validation. Finally, we show some IDF curves modeled with the best model.

## 3.1. Assessment of EGPD Goodness of Fit

The first issue is to investigate whether EGPD is an appropriate model for the precipitation data at hand. To check this, we fitted the model at each station and for each duration, independently. We call this EGPD model fitted on each data separately as the “base” model. We then assess the quality of the resulting fits by computing the NRMSE given in Equation 20. The seasonal boxplots of the score for each duration are shown in Figure 5. The higher the score, the better the model.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Boxplots of $1 - \text{NRMSE}_\epsilon(d)$ versus duration for the base extended generalized Pareto distribution model, that is, fitted on data of each duration separately. Each boxplot contains 81 points, with each point corresponding to one station.</figcaption></figure>

In spring and summer, the quality of the fit is less good for durations lower than 2 hr. In winter, on the other hand, the fit is less good for $d = 48$ and 72 hr. Overall, more than 74% of the scores fall above 0.9 and 96% above 0.8. We, therefore, consider the EGPD to be a reasonable model for the data.

The fitted shape parameter $\xi$ with respect to duration is shown on Figure 6. Each boxplot contains 81 values, one for each station. We can observe strong dependence of this parameter on duration, especially in summer. For this season, while 75% of the stations have a $\xi > 0.17$ for $d = 1$ hr, only 25% have $\xi > 0.06$ at $d = 24$ hr. In winter, however, the dependence is not very strong, as judged by the large variability of the boxplots.

## 3.2. Comparison of Models

Results of the model comparison are presented under two frameworks, first in calibration, and then second in cross-validation based on split sampling.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Boxplots of the fitted $\xi$ versus duration obtained with the base model. Each boxplot contains 81 points, with each point corresponding to one station.</figcaption></figure>

### 3.2.1. Evaluation in Calibration

#### 3.2.1.1. Normalized Root Mean Square Error

Figure 7 presents the seasonal boxplots of the 1-NRMSE for the 10 IDF models and the base model. Each boxplot contains 1,053 points, summarizing the score over 81 stations and 13 durations. In the case of the base model (in yellow), the scores are the same as those in Figure 5, but here we merge the scores for all the durations together.

For all seasons, the two data-driven IDF models, $IDF_{DD_{local}}$ and $IDF_{DD_{global}}$, always show the best performance compared to the others. When looking at the two, the $IDF_{DD_{global}}$ generally outperforms the $IDF_{DD_{local}}$. This means that the global fitting of the model improves the estimation performance compared to the simple interpolation of the locally estimated parameters.

Comparing the $IDF_{ss}$ and the $IDF_{kouts}$ (white vs. red boxplots), the results show that for all seasons, the $IDF_{kouts}$ has a better performance compared to the $IDF_{ss}$. Recall that the two models differ by the additional parameter $\theta$ in the former to account for curvature for short durations.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Boxplots of the (1-NRMSE) in calibration. Each boxplot contains 1,053 points, each point corresponding to one station and duration.</figcaption></figure>

Allowing for $\xi = f(d)$ (models with subscript _$\xi(d)$) increased the performance of the models mainly in summer, where all the models without this addition showed very poor performance. For the other seasons, the gain in performance is not as pronounced.

Finally, the models allowing for scaling break (those with subscript _TR), show improved performance compared to those with the single regime for all the seasons, except summer (e.g., IDF$_{ss}$ vs. IDF$_{ss\_TR}$, i.e., the white and violet boxplots).

We note here that we used the NRMSE (Figure 7) to measure the accuracy over all non-zero intensities. This is because we are using the EGPD model, which is supposed to model correctly all the non-zero intensities. However, we also computed the NRMSE on extremes only, which we define as the exceedances over the 98% quantile overall intensities (including zeros). The normalization was done by the average of the exceedances. The result (not shown) maintains the same performance ranking order of the models as in Figure 7.

### 3.2.1.2. Akaike Information Criteria

The result in Figure 7 as measured by the NRMSE assesses the in-sample accuracy of the models but does not reward parsimony in terms of the number of model parameters. A natural question to ask is whether the additional
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Akaike Information Criteria (AIC) results for the four seasons. Each panel shows the percentage of stations (y-axis) over which a model is ranked (x-axis) according to the AIC criteria. The ranks are such that “1” is the best and “10” is the worst.</figcaption></figure>

performance is worth the additional complexity. To answer this, we compute the AIC for each model at each of the 81 stations. Each time, we rank the models from the best (rank = 1, smallest AIC) to the least (rank = 10, largest AIC). Figure 8 shows results in the form of stacked-bar chars for the four seasons. The horizontal axis shows the rank from 1 to 10, and the vertical axis shows the percentage of stations over which a model is ranked. For instance, in winter, IDF<sub>DDglobal</sub> is the best (rank = 1) in 82.7% of the stations, while the base model is the best in 17.3% of the stations.

The result shows that for all the seasons, the first and second rank is almost exclusively shared between the IDF<sub>DDglobal</sub> and the base model. To summarize the results, Table 2 gives the ranking of the models for the four seasons. The model with the highest percentage is selected for each rank (most likely model at each rank). For all the seasons, the Data-driven (IDF<sub>DDglobal</sub>) is the best, the base model following behind. The results also revealed the relevance of the models allowing for shape parameters to vary with duration only in summer. For the other three seasons, however, the gain in performance is not worth the additional parameter modeling of the shape parameter as a function of duration. Comparing the models with a constant shape parameter, while the IDF<sub>kouts</sub> is always the most parsimonious, the simple-scaling (IDF<sub>ss</sub>) remained the worst.
<table><caption>**Table 2**
Summary of the AIC Ranking of the Models in the Four Seasons</caption><thead><tr><th>Model</th><th>No. of parameters</th><th>Winter</th><th>Spring</th><th>Summer</th><th>Autumn</th></tr></thead><tbody><tr><td>Base</td><td>39</td><td>2</td><td>2</td><td>2</td><td>2</td></tr><tr><td>IDF<sub>DDglobal</sub></td><td>10</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>IDF<sub>ss</sub></td><td>4</td><td>6</td><td>6</td><td>10</td><td>6</td></tr><tr><td>IDF<sub>ss_ξ(d)</sub></td><td>5</td><td>10</td><td>10</td><td>6</td><td>10</td></tr><tr><td>IDF<sub>kouts</sub></td><td>5</td><td>3</td><td>3</td><td>7</td><td>3</td></tr><tr><td>IDF<sub>kouts_ξ(d)</sub></td><td>6</td><td>9</td><td>7</td><td>3</td><td>9</td></tr><tr><td>IDF<sub>ss_TR</sub></td><td>6</td><td>5</td><td>4</td><td>8</td><td>5</td></tr><tr><td>IDF<sub>ss_TR_ξ(d)</sub></td><td>7</td><td>8</td><td>9</td><td>4</td><td>7</td></tr><tr><td>IDF<sub>kouts_TR</sub></td><td>7</td><td>4</td><td>5</td><td>9</td><td>4</td></tr><tr><td>IDF<sub>kouts_TR_ξ(d)</sub></td><td>8</td><td>7</td><td>8</td><td>5</td><td>8</td></tr></tbody></table>

*Note.* The model with rank 1 is the best, that is, it has the smallest AIC for most of the stations in that season. A rank of 10 indicates the worst-performing model.

In summary, the *IDF<sub>DDglobal</sub>* which had the best in-sample performance among the IDF models (as measured by the NRMSE in Figure 7), remained the best across all the seasons based on the AIC criteria. Although our focus has been to compare the IDF models that consider the linkage of parameters with duration, we still included the base model in the comparison of the AIC. The result here showed that the base model (with 39 parameters) is less parsimonious compared to the *IDF<sub>DDglobal</sub>* for most of the stations.

Finally, it was not possible to consider the *IDF<sub>DDlocal</sub>* here, because, unlike the 10 models that are based on maximum likelihood estimation, this model involves a two-step estimation process. The first is a maximum likelihood, followed by a second step to estimate the link between the parameters and duration through least squares and segmented regressions. Nevertheless, this is not an issue because this model is based on the same principle as the *IDF<sub>DDglobal</sub>* and simply differs by the estimation method. Thus, comparing the two *IDF<sub>DDglobal</sub>* versus *IDF<sub>DDlocal</sub>* can be seen as comparing two estimation methods within the same model, and according to the NRMSE results in Figure 7, better performance is achieved with the *IDF<sub>DDglobal</sub>* compared to the *IDF<sub>DDlocal</sub>*.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Boxplots of the regional cross-validation criteria, NRMSE_CV, FF, and SPAN. For the first two criteria, each boxplot contains 2 × 130 points, corresponding to one regional score for each of the 13 durations and 10 repetitions of the split sampling. For the SPAN, each boxplot contains 130 points. The optimal value for each criterion is equal to 1.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Simple-scaling (IDF<sub>ss</sub>) and data-driven (IDF<sub>DDglobal</sub>) curves (a) in summer at a station in Zurich (north-east). (b) In autumn at a station in Locarno (Ticino area in the south). The curves are for the return periods *T* = 2, 5, 10, 40, and 100 years. The points are the empirical quantiles corresponding to *T* = 2, 5, and 10 years. The envelopes represent the 95% confidence bounds obtained by block bootstrap.</figcaption></figure>

### 3.2.2. Evaluation in Cross-Validation

The split-sampling procedure allows for the comparison of the models in a cross-validation framework. We use three regional criteria: NRMSE_CV, FF, and SPAN (see Section 3.5.2), to enable the comparison of the models based on their predictive capabilities. We want to select a model, which in addition to being able to fit the data used to train it, is able to perform reliably and robustly in the presence of new data.
In the following, we present the results in three paragraphs, first according to the reliability/accuracy of the model in predicting all the observations as measured by NRMSE_CV, then the reliability in predicting the maxima as measured by the FF criterion, and finally, the robustness of the model in predicting the 100-year return level as measured by SPAN100. Figure 9 presents the results for the four seasons. For all the criteria, the model with a regional score equal to 1 is the best model.

For all seasons, the NRMSE_CV shows the data-driven models, specifically the *IDF<sub>DDglobal</sub>*, to be the most accurate/reliable in predicting the entire observations compared to the other models. In winter, however, the difference in the performance of the models is not very clear. Looking at the summer results, the models without accounting for $\xi = f(d)$ always have the worst performance.

In terms of the FF criterion, the best performance in predicting the maxima in winter is shown by the *IDF<sub>ss_TR</sub>* model. In fact, all the models with no allowance for $\xi = f(d)$ happen to be the most reliable models in this season. The converse is however true in the case of the remaining seasons. In summer, while *IDF<sub>DDlocal</sub>* is the best model, *IDF<sub>ss_ξ(d)</sub>* is the best in spring and autumn.

The robustness criteria, SPAN100 shows the models with no allowance for $\xi = f(d)$ to be the most robust models. An exception to this is in summer, where the *IDF<sub>DDglobal</sub>* model is the most robust model. Also, higher robustness is found for the models not accounting for $\xi = f(d)$ compared to their counterparts, for example, *IDF<sub>ss</sub>* versus *IDF<sub>ss_ξ(d)</sub>*. This is despite the fact that the former performs poorly in calibration, and is the least performing according to the other cross-validation criteria of reliability. This confirms the previous comments of Garavaglia et al. (2011) that a robust model can completely fail to model/predict the data. Hence the robustness criteria should only be used alongside other reliability criteria, such that the most robust model is only selected among models of similar reliability.

To summarize the results, the best IDF model should perform well in calibration, and should not be very sensitive to the data used to train it. In calibration, the data-driven model *IDF<sub>DDglobal</sub>* showed the best performance compared to all the other nine models, it also remains accurate and reliable at predicting the entire observations in the split-sampling cross-validation (as measured by the NRMSE_CV), especially in summer. This is an important feature since we are interested in the complete range of intensities. Finally, it generally showed more robustness compared to the other models of similar reliability.

### 3.3. IDF Curves

Figure 10a shows the IDF curves from two models, *IDF<sub>ss</sub>* and *IDF<sub>DDglobal</sub>* along with their 95% CI in summer, at a station in Zurich which is located in the Northeast of Switzerland. In this region, summer is the main season of heavy rainfall. As a reminder, the *IDF<sub>ss</sub>* allows scaling only in the scale parameter, $\sigma$ of the EGPD, the other two parameters ($\kappa$ and $\xi$), are independent of duration. The *IDF<sub>DDglobal</sub>* on the other hand, allows each of the three parameters to vary with duration. The curves are for return periods $T = 2, 5, 10, 50$, and 100 years, while points are the empirical levels for $T = 2, 5$, and 10 years. The *IDF<sub>ss</sub>* performed poorly at predicting the empirical quantiles. The curves modeled by the *IDF<sub>DDglobal</sub>* on the other hand, are in agreement with the empirical levels. Similar IDF curves for autumn are shown in Figure 10b for a station in Locarno which is located in the Ticino area in the south of Switzerland. The Ticino area is subject to the heaviest precipitation compared to the other regions in Switzerland. Again, the *IDF<sub>DDglobal</sub>* is able to model the empirical levels correctly for both the short and long durations.

The 95% confidence bounds for the return periods $T = 2, 5, 10, 50$, and 100 years in both stations are shown in Figure 11 along with the empirical levels for $T = 2, 5$, and 10 years (broken lines). For comparison sake, we also include those of the base model, that is, where EGPD is fitted to the data of each duration separately (without linkage of parameters to duration). It can be seen that narrower CIs are obtained with the simple scaling (*IDF<sub>ss</sub>*), and data-driven *IDF<sub>DDglobal</sub>* compared to the base model. This is because the base model uses less amount of data in its inference compared to the other two models where all the data of the 13 durations are pooled together in the estimation. It is expected that the width of the CI would be narrower as the available data for estimation increases. The *IDF<sub>ss</sub>* has six parameters less than those of the *IDF<sub>DDglobal</sub>*, so a natural question is whether the CI will be wider in the latter model. Indeed for levels where both models predicted the empirical level correctly ($d = 1$ hr in the case of Locarno, second row of Figure 11), the CI is a bit narrower in the case of the *IDF<sub>ss</sub>*. But looking at
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Return level estimates (points) and their 95% confidence interval bounds (lines) for return levels *T* = 2, 5, 10, 50, and 100 years estimated with the simple-scaling (IDF<sub>ss</sub>), data-driven (IDF<sub>DDglobal</sub>) and the base model (no linkage of extended generalized Pareto distribution parameters with duration). The confidence bounds were obtained using block bootstrap. The top row is for a station in Zurich (northeast) in summer for *d* = 1 and 24 hr. The second row is for a station in Locarno (Ticino area in the south) in autumn for *d* = 1 and 24 hr. In all cases, the broken lines in red, black, and yellow are the empirical levels for *T* = 2, 5, and 10 years respectively.</figcaption></figure>

the same station for *d* = 24 hr, the bounds are wider, in addition to the bias in predicting the empirical level. Of course, a narrower CI is only relevant if it contains the empirical level.

In Figure 10, the curves of the simple-scaling model (IDF<sub>ss</sub>) are not parallel. This behavior resulted from the defi-
nition of IDF models for non-zero precipitation in Equation 5. From this equation, we see that the *T*-year return
level is defined as ($1 - \frac{1}{T \times \delta_d}$). The term $\delta_d$, representing the average number of non-zero precipitation per year
varies across the durations leading to a non-constant slope for the different curves.

We finally show, in Figures 12 and 13, respectively, the seasonal 100-year return level maps for *d* = 1 and 24 hr. The levels in the maps can be interpreted as the levels which are expected to be exceeded every 100 seasons. The levels were obtained with the best performing model, that is, IDF<sub>DDglobal</sub>. Looking at the return levels for *d* = 1 hr (Figure 12), we see that the levels in winter are the lowest, with no specific spatial pattern or variability. In spring, the levels in the north and Ticino starts to increase. Summer has the highest levels, and similar levels are obtained all along the plateau in the north. In autumn, while the levels in the north are comparable to those in
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Map of Switzerland showing the seasonal 100-year return level in mm/hr for *d* = 1 hr. Levels predicted with IDF<sub>DDglobal</sub>.</figcaption></figure>

winter, those in Ticino are comparable to those in summer. A different spatial pattern is however observed for the 100-year return level for *d* = 24 hr. Specifically for summer, the levels in the plateau are lower than those along the northern alpine rim (Prealps). The exhibited spatial pattern of the levels produced by this model is similar to those observed in earlier studies (see Fukutome et al. (2015) for the hourly, and Haruna et al. (2022) for the 24 hr precipitation).

**3.4. Discussion**

In the following paragraphs, we briefly discuss some of our choices in terms of the functional forms of the data-driven models, taking into account the varying shape parameter with duration, and the issue of scaling break in the data.

First, for the data-driven models, we limited our choice of functional relationships to simple parametric models, specifically to piece-wise linear models. Other choices would be possible such as smooth regression splines (e.g., Youngman, 2019, 2020). This choice has its advantage and drawback. The advantage is that the splines are able to automatically adjust to fit any form of relationship. The main drawback is that it is inherently non-parametric, and so the mapping of the IDF models, to allow predictions at ungauged locations, is not directly possible. One can only map the three EGPD parameters for a particular duration. For instance, for 13 durations, this means 3 × 13 = 39 maps. For our choice of linear functions, 10 parameters are able to describe the IDF curves at each station, and hence 10 maps for the whole area under study.

Regarding the variation of the shape parameter ξ with respect to duration *d*, some earlier studies did observe or discuss it (e.g., Ulrich et al., 2021; Veneziano et al., 2007). They however did not model it, either due to the weak form of the relationship, or because the IDF model did not allow for it. Here, especially in summer, we found very strong dependence, and the results have shown that taking it into account is invaluable. For the other seasons, however, the AIC results showed that it is more parsimonious to let the parameter be independent of duration.
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Map of Switzerland showing the seasonal 100-year return level in mm/hr for *d* = 24 hr. Levels predicted with IDF<sub>DDglobal</sub>.</figcaption></figure>

We note that the strength of the variation of the shape parameter with duration is different in the four seasons. In our case, it is strongest in summer and weakest in winter. A possible explanation is that the underlying precipitation formation mechanism responsible for the formation of short and long-duration precipitation is different in the four seasons. For instance, in winter, the same frontal system is responsible for both short and long-duration events, and hence the variation of the shape parameter with duration is weak. In summer where the linkage is strong, different systems are responsible for the short-duration intensities (e.g., due to short convective events) and long-lasting duration intensities (e.g., due to frontal systems). The other seasons (spring and autumn) present a mixture of the two behaviors.

Finally, we observed a break in the scaling relationship of the EGPD scale parameter. The existence of scaling regimes has been investigated by Fraedrich and Larnder (1993) and they linked it to atmospheric processes such as the structure of frontal systems, the diurnal cycle, the seasonal periodicity, and climatic variability. Hence, break points in the regime observed at stations could be explained by the possible transition from one rainfall system to another. For example, from short convective events at short duration to frontal systems at longer duration. For breaks at very short durations in the range of around 1 hr, however, it is generally attributed to being artificial, imposed by the resolution of the gauge (measurement precision errors). Despite this, the identification of scaling breaks involves some level of uncertainty depending on the method used. For instance, Paschalis (2013) tried to identify scaling breaks in Switzerland using the same 10 min data sets as in our case. He used spectral density analysis based on the Fourier transform, and another based on wavelet decomposition. While both methods agree on the existence of scaling breaks, the exact times of the breaks were different in both. The author also mentioned that no clear seasonal or regional pattern was observed. In our case, we rely on the estimated EGPD scale parameter to identify the scaling breaks. Other authors relied on moment analysis (plot of moment vs. duration in log-log scale for different orders) to identify breaks in scaling (Bougadis & Adamowski, 2006; Nhat et al., 2008). In reference to this, we do not rule out the possible effect of the EGPD model on the observed scaling breaks. A detailed study of its own might be required to study and characterize/interpret the scaling breaks in the study area since any detailed interpretation might be premature.
4. Conclusions

Our aim in this paper was to build IDF curves using all the non-zero precipitation data in Switzerland. To achieve this, we used the EGPD model as the parametric model for the precipitation intensities. The literature presents various approaches to link the different durations together in IDF curves. We considered three of these approaches to build the IDF curves while using the EGPD as the parametric model. The first approach is based on the scale invariance theory, where IDF curves are built based on the scaling behavior of precipitation. The second approach is based on the general IDF formulation of Koutsoyiannis et al. (1998), which generalizes the various traditional IDF formulations. The last approach is called the data-driven approach, where each parameter can vary with duration, and the form of the relationship is empirically determined by the data at hand.

We started from these three approaches and added some extensions to account for scaling break and varying
shape parameter with duration. We ended up with a total of 10 IDF models. We then compared them, first in
calibration, and then in a split-sample cross-validation approach.

The results showed that, given the EGPD as the parametric model, the data-driven IDF-EGPD, particularly the IDF<sub>DDglobal</sub>, is the best model for the data at hand. The AIC also showed this model to be more parsimonious compared to the base model where no linkage of parameter with duration is considered. The IDF curves based on simple-scaling and the general formulation of Koutsoyiannis et al. (1998), did not perform as efficiently even with the added extensions in terms of scaling break and in the way the shape parameter varies with duration. The fact that the simple-scaling IDF models performed poorly in summer confirms the previous findings of Molnar and Burlando (2008) and Paschalis (2013) that in Switzerland, precipitation in summer shows multiscaling behavior.

Although our work focused on Switzerland, the data-driven approach can be applied everywhere, especially for regions where high-resolution data are available. This is because the approach considers the empirical relationship between the model parameters and duration, and is not constrained by the hypothesis of scaling across durations. Our result also showed that it is possible to model the linkage of the shape parameter with duration in IDF curve modeling. Moreover, the study highlights the need to explore in detail the empirical relationship of model parameters with duration before the application of any of the widely used IDF construction approaches. Finally, we showed that it is possible to use the EGPD for IDF curve modeling.

In terms of perspectives, it would be interesting to produce maps of the parameters to allow for predictions at ungauged sites. This could be achieved by simple interpolation of the local IDF parameters as done by Blanchet et al. (2016), or through quantile regression methods (Ouali & Cannon, 2018), or by global estimation using spatial covariates (e.g., Ulrich et al., 2020). Another possibility is to use a regionalization technique, such as the method of Hosking and Wallis (2005) and then interpolate the index flood to allow predictions at the ungauged sites (e.g., Mascaro, 2020).

It should be mentioned that throughout this work, we estimated the IDF models through the independence like-
lihood, thus omitting the correlation between different timesteps and durations. In the case of annual maxi-
mum series (GEV), Nadarajah et al. (1998) have modeled the dependence between the different time steps using
multivariate extreme value distributions, and Tyralis and Langousis (2019) followed suit by using max-stable
processes. Later, Jurado et al. (2020) investigated the impact of accounting for this dependence in extremes and
showed that there is little gain in performance, in addition to the added complexity of using max-stable processes.
Since all these authors focused on the distributions that are based on extreme data only, an interesting perspective
will be to investigate the effect of accounting for the dependence in IDF curve modeling for the case of the EGPD
that uses all the non-zero data.

Finally, consideration of the effect of climate change in building IDF curves is invaluable. For instance, Cheng and AghaKouchak (2014) showed that by neglecting nonstationarity in modeling IDF curves, there could be up to 60% underestimation of extreme precipitation, especially for short durations. It would therefore be interesting to model the curves while accounting for a warmer climate (e.g., Cheng & AghaKouchak, 2014; Kristvik et al., 2019; Mirhosseini et al., 2013; Ouarda et al., 2019; Ragno et al., 2018).

Appendix A: Variables and Their Meaning

Table A1.
<table><caption>**Table A1**
List of Variables and Their Meaning</caption><thead><tr><th></th><th>Variable</th><th>Meaning</th><th>Range</th><th>Unit</th></tr></thead><tbody><tr><td>1</td><td>H</td><td>Power-law scaling exponent</td><td>0 &lt; H &lt; 1</td><td>–</td></tr><tr><td>2</td><td>θ</td><td>Duration offset</td><td>θ &gt; 0</td><td>hr</td></tr><tr><td>3</td><td>d<sub>0</sub></td><td>Reference duration</td><td>d<sub>0</sub> = 1 in Simple scaling d<sub>0</sub> = 1 – θ in General IDF formulation</td><td>hr</td></tr><tr><td>4</td><td>κ<sub>d0</sub></td><td>EGPD kappa parameter for d<sub>0</sub></td><td>κ<sub>d0</sub> &gt; 0</td><td>–</td></tr><tr><td>5</td><td>σ<sub>d0</sub></td><td>EGPD scale parameter for d<sub>0</sub></td><td>σ<sub>d0</sub> &gt; 0</td><td>mm/hr</td></tr><tr><td>6</td><td>ξ<sub>d0</sub></td><td>EGPD shape parameter for d<sub>0</sub></td><td>ξ<sub>d0</sub> ≥ 0</td><td>–</td></tr><tr><td>7</td><td>K, K*</td><td>Duration of scaling break</td><td>K &gt; 0, K* &gt; 0</td><td>hr</td></tr><tr><td>8</td><td>H<sub>1</sub>, H<sub>2</sub></td><td>1st and 2nd slopes in two-scaling regimes</td><td>0 &lt; H<sub>1</sub> &lt; 1, 0 &lt; H<sub>2</sub> &lt; 1</td><td>–</td></tr></tbody></table>

## Appendix B: EGPD Parameters Versus Duration

### B1. Parametric Models for the EGPD Shape Parameter Versus Duration

* constant:

$$ \xi_d = \frac{1}{n} \sum_{j=1}^{n} \xi_{d,j} \qquad (B1) $$

* linear:

$$ \xi_d = a_{\xi} + b_{\xi} d \qquad (B2) $$

* log-linear:

$$ \log(\xi_d) = a_{\xi} + b_{\xi} d \qquad (B3) $$

* linear-log:

$$ \xi_d = a_{\xi} + b_{\xi} \log(d) \qquad (B4) $$

* log-log:

$$ \log(\xi_d) = a_{\xi} + b_{\xi} \log(d) \qquad (B5) $$

### B2. Plots of the EGPD Parameters Versus Duration

In this section, we explore the relationship between the three EGPD parameters with respect to duration. Two stations located at Zurich (KLO) and Robbia in Graubünden (ROB) are given as examples. For each season, the fitted parameter from the base model is plotted against the duration *d*.
<figure><img src="image_17.png" /><figcaption>**Figure B1.** Relationship between the extended generalized Pareto distribution (EGPD) $\kappa$ parameter and duration $d$ for a station in Zurich (KLO). To obtain this, we fitted EGPD to the data of each of the 13 durations separately by maximum likelihood. The black points are the fitted $\kappa$. The broken lines in red represent the 95% confidence interval for the fitted model (log-log TR) see Equation 11.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure B2.** Same as in Figure B1, but for the extended generalized Pareto distribution scale parameter $\sigma$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure B3.** Same as in Figure B1, but for the extended generalized Pareto distribution shape parameter $\xi$. The equations for the fitted models are given in Appendix B. The broken lines represent the 95% confidence interval for the linear-log model.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure B4. Same as in Figure B1, but for a station at Robbia in Graubünder (ROB).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure B5. Same as in Figure B4, but for the extended generalized Pareto distribution scale parameter σ.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure B6.** Same as in Figure B4, but for the extended generalized Pareto distribution shape parameter ξ. The equations for the fitted models are given in Appendix B. The broken lines represent the 95% confidence interval for the linear-log model.</figcaption></figure>

## Data Availability Statement

The data used in this study are maintained by the Swiss Federal Office of Meteorology and Climatology, MeteoSwiss, and can be freely obtained from IDAWEB (2021). The software for this research is freely available as an R package (Haruna, 2023) with GPL-3 license.

## References

Agbazo, M. N., Koto N'Gobi, G., Kounouhewa, B., Alamou, E., Afouda, A., & Akpo, A. (2016). Estimation of IDF curves of extreme rainfall by simple scaling in Northern Oueme Valley, Benin Republic (West Africa). *Earth Sciences Research Journal*, 20(1), 1-7. https://doi.org/10.15446/esrj.v20n1.49405

Akaike, H. (1974). A new look at the statistical model identification. *IEEE Transactions on Automatic Control*, 19(6), 716-723. https://doi.org/10.1109/TAC.1974.1100705

Ben-Zvi, A. (2009). Rainfall intensity-duration-frequency relationships derived from large partial duration series. *Journal of Hydrology*, 367(1-2), 104-114. https://doi.org/10.1016/j.jhydrol.2009.01.007

Bernard, M. M. (1932). Formulas for rainfall intensities of long duration. *Transactions of the American Society of Civil Engineers*, 96, 592-606. https://doi.org/10.1061/taceat.0004323

Blanchet, J., Ceresetti, D., Molinié, G., & Creutin, J.-D. (2016). A regional GEV scale-invariant framework for Intensity-Duration-Frequency analysis. *Journal of Hydrology*, 540, 82-95. https://doi.org/10.1016/j.jhydrol.2016.06.007

Blanchet, J., Paquet, E., Vaittinaida Ayar, P., & Penot, D. (2019). Mapping rainfall hazard based on rain gauge data: An objective cross-validation framework for model selection. *Hydrology and Earth System Sciences*, 23(2), 829-849. https://doi.org/10.5194/hess-23-829-2019

Blanchet, J., Touati, J., Lawrence, D., Garavaglia, F., & Paquet, E. (2015). Evaluation of a compound distribution based on weather pattern subsampling for extreme rainfall in Norway. *Natural Hazards and Earth System Sciences*, 15(12), 2653-2667. https://doi.org/10.5194/nhess-15-2653-2015

Bougadis, J., & Adamowski, K. (2006). Scaling model of a rainfall intensity-duration-frequency relationship. *Hydrological Processes*, 20(17), 3747-3757. https://doi.org/10.1002/hyp.6386

Burlando, P., & Rosso, R. (1996). Scaling and multiscaling models of depth-duration-frequency curves for storm precipitation. *Journal of Hydrology*, 187(1), 45-64. https://doi.org/10.1016/S0022-1694(96)03086-7

Carreteras, M. T. (1987). *Cálculo hidrometeorológico de caudales máximos en pequeñas cuencas naturales*. Textos de la Dirección General de Carreteras, (12), 1987.

## Acknowledgments

This research has been supported by the Bundesamt für Energie (Grant SI/502150-01) and the Bundesamt für Umwelt (Grant SI/502150-01), through the *Extreme floods in Switzerland* project.
Chang, K. B., Lai, S. H., & Othman, F. (2016). Comparison of annual maximum and partial duration series for derivation of rainfall intensity-duration-frequency relationships in Peninsular Malaysia. *Journal of Hydrologic Engineering*, 21(1), 05015013. https://doi.org/10.1061/(ASCE)HE.1943-5584.0001262

Cheng, L., & AghaKouchak, A. (2014). Nonstationary precipitation intensity-duration-frequency curves for infrastructure design in a changing climate. *Scientific Reports*, 4(1), 1–6. https://doi.org/10.1038/srep07093

Chow, V. (1962). *Hydrologic determination of waterway areas for drainage structures in small drainage basins*, Engrg (p. 11). Experimental Station, University of Illinois.

Coles, S. (2001). *An introduction to statistical modeling of extreme values*. Springer.

Courty, L. G., Wilby, R. L., Hillier, J. K., & Slater, L. J. (2019). Intensity-duration-frequency curves at the global scale. *Environmental Research Letters*, 14(8), 084045. https://doi.org/10.1088/1748-9326/ab370a

Davison, A. C., Padoan, S. A., & Ribatet, M. (2012). Statistical modeling of spatial extremes. *Statistical Science*, 27(2), 161–186. https://doi.org/10.1214/11-STS376

Eicher, C., & Krejci, V. (1996). A new rainfall data system for urban hydrology in Switzerland. *Atmospheric Research*, 42(1), 177–198. https://doi.org/10.1016/0169-8095(95)00062-3

Evin, G., Blanchet, J., Paquet, E., Garavaglia, F., & Penot, D. (2016). A regional model for extreme rainfall based on weather patterns subsampling. *Journal of Hydrology*, 541, 1185–1198. https://doi.org/10.1016/j.jhydrol.2016.08.024

Evin, G., Favre, A.-C., & Hingray, B. (2018). Stochastic generation of multi-site daily precipitation focusing on extreme events. *Hydrology and Earth System Sciences*, 21(1), 655–672. https://doi.org/10.5194/hess-22-655-2018

Fauer, F. S., Ulrich, J., Jurado, O. E., & Rust, H. W. (2021). Flexible and consistent quantile estimation for intensity-duration-frequency curves. *Hydrology and Earth System Sciences*, 25(12), 6479–6494. https://doi.org/10.5194/hess-25-6479-2021

Fraedrich, K., & Larnder, C. (1993). Scaling regimes of composite rainfall time series. *Tellus*, 45(4), 289–298. https://doi.org/10.1034/j.1600-0870.1993.t01-3-00004.x

Frei, C., & Schär, C. (1998). A precipitation climatology of the Alps from high-resolution rain-gauge observations. *International Journal of Climatology: A Journal of the Royal Meteorological Society*, 18(8), 873–900. https://doi.org/10.1002/(sici)1097-0088(19980630)18:8<873::aid-joc255>3.0.co;2-9

Froidevaux, P., Schwanbeck, J., Weingartner, R., Chevalier, C., & Martius, O. (2015). Flood triggering in Switzerland: The role of daily to monthly preceding precipitation. *Hydrology and Earth System Sciences*, 19(9), 3903–3924. https://doi.org/10.5194/hess-19-3903-2015

Fukutome, S., Liniger, M. A., & Süveges, M. (2015). Automatic threshold and run parameter selection: A climatology for extreme hourly precipitation in Switzerland. *Theoretical and Applied Climatology*, 120(3–4), 403–416. https://doi.org/10.1007/s00704-014-1180-5

Garavaglia, F., Lang, M., Paquet, E., Gailhard, J., Garçon, R., & Renard, B. (2011). Reliability and robustness of rainfall compound distribution model based on weather pattern sub-sampling. *Hydrology and Earth System Sciences*, 15(2), 519–532. https://doi.org/10.5194/hess-15-519-2011

Ghanmi, H., Bargaoui, Z., & Mallet, C. (2016). Estimation of intensity-duration-frequency relationships according to the property of scale invariance and regionalization analysis in a mediterranean coastal area. *Journal of Hydrology*, 541, 38–49. https://doi.org/10.1016/j.jhydrol.2016.07.002

Giannakaki, P., & Martius, O. (2015). Synoptic-scale flow structures associated with extreme precipitation events in northern Switzerland. *International Journal of Climatology*, 40(6), 19–2515. https://doi.org/10.1002/joc.4508

Gupta, V. K., & Waymire, E. (1990). Multiscaling properties of spatial rainfall and river flow distributions. *Journal of Geophysical Research: Atmospheres*, 95(D3), 1999–2009. https://doi.org/10.1029/JD095iD03p01999

Gupta, V. K., & Waymire, E. C. (1993). A statistical analysis of mesoscale rainfall as a random cascade. *Journal of Applied Meteorology and Climatology*, 32(2), 251–267. https://doi.org/10.1175/1520-0450(1993)032(0251:ASAOMR)2.0.CO;2

Harris, D., Seed, A., Menabde, M., & Austin, G. (1997). Factors affecting multscaling analysis of rainfall time series. *Nonlinear Processes in Geophysics*, 4(3), 137–156. https://doi.org/10.5194/npg-4-137-1997

Haruna, A. (2023). egpdidf: Modeling intensity-duration-frequency (idf) curves using the extended generalized pareto distribution (egpd) [software]. Zenodo. https://doi.org/10.5281/zenodo.7828750

Haruna, A., Blanchet, J., & Favre, A.-C. (2022). Performance-based comparison of regionalization methods to improve the at-site estimates of daily precipitation. *Hydrology and Earth System Sciences*, 26(10), 2797–2811. https://doi.org/10.5194/hess-26-2797-2022

Hosking, J. R. M., & Wallis, J. R. (2005). *Regional frequency analysis: An approach based on L-moments*. Cambridge University Press.

IDAWEB. (2021). The data portal of meteoswiss for research and teaching [Dataset]. Retrieved 2021-12-10, from https://gate.meteoswiss.ch/idaweb/login.do

Innocenti, S., Mailhot, A., & Frigon, A. (2017). Simple scaling of extreme precipitation in North America. *Hydrology and Earth System Sciences*, 21(11), 5823–5846. https://doi.org/10.5194/hess-21-5823-2017

Isotta, F. A., Frei, C., Weilguni, V., Perčec Tadić, M., Lassègues, P., Rudolf, B., et al. (2014). The climate of daily precipitation in the alps: Development and analysis of a high-resolution grid dataset from pan-alpine rain-gauge data: Climate of daily precipitation in the alps. *International Journal of Climatology*, 34(5), 1657–1675. https://doi.org/10.1002/joc.3794

Jurado, O. E., Ulrich, J., Scheibel, M., & Rust, H. W. (2020). Evaluating the performance of a max-stable process for estimating intensity-duration-frequency curves. *Water*, 12(12), 3314. https://doi.org/10.3390/w12123314

Katz, R. W., Parlange, M. B., & Naveau, P. (2002). Statistics of extremes in hydrology. *Advances in Water Resources*, 25(8–12), 1287–1304. https://doi.org/10.1016/S0309-1708(02)00056-8

Koutsoyiannis, D., Kozonis, D., & Manetas, A. (1998). A mathematical framework for studying rainfall intensity-duration-frequency relationships. *Journal of Hydrology*, 206(1–2), 118–135. https://doi.org/10.1016/S0022-1694(98)00097-3

Kristvik, E., Johannessen, B. G., & Muthanna, T. M. (2019). Temporal downscaling of IDF curves applied to future performance of local storm-water measures. *Sustainability*, 11(5), 1231. https://doi.org/10.3390/su11051231

Langousis, A., Mamalakis, A., Puliga, M., & Deidda, R. (2016). Threshold detection for the generalized Pareto distribution: Review of representative methods and application to the NOAA NCDC daily rainfall database. *Water Resources Research*, 52(4), 2659–2681. https://doi.org/10.1002/2015WR018502

Langousis, A., & Veneziano, D. (2007). Intensity-duration-frequency curves from scaling representations of rainfall. *Water Resources Research*, 43(2), W02422. https://doi.org/10.1029/2006WR005245

Le Gall, P., Favre, A.-C., Naveau, P., & Prieur, C. (2022). Improved regional frequency analysis of rainfall data. *Weather and Climate Extremes*, 36, 100456. https://doi.org/10.1016/j.wace.2022.100456

Lima, M. I. P. d. (1998). *Multifractals and the temporal structure of rainfall*. PhD Thesis. Wageningen Agricultural University.
Madsen, H., Rosbjerg, D., & Harremoëes, P. (1995). Application of the Bayesian approach in regional analysis of extreme rainfalls. *Stochastic Hydrology and Hydraulics*, 9(1), 77–88. https://doi.org/10.1007/BF01581759

Mascaro, G. (2020). Comparison of local, regional, and scaling models for rainfall intensity-duration-frequency analysis. *Journal of Applied Meteorology and Climatology*, 59(9), 1519–1536. https://doi.org/10.1175/jamc-d-20-0094.1

Mélése, V., Blanchet, J., & Molinié, G. (2018). Uncertainty estimation of intensity-duration-frequency relationships: A regional analysis. *Journal of Hydrology*, 558, 579–591. https://doi.org/10.1016/j.jhydrol.2017.07.054

Menabde, M., Seed, A., & Pegram, G. (1999). A simple scaling model for extreme rainfall. *Water Resources Research*, 35(1), 335–339. https://doi.org/10.1029/1998WR900012

Meylan, P., Favre, A.-C., & Musy, A. (2012). *Predictive hydrology: A frequency analysis approach*. CRC Press.

Mirhosseini, G., Srivastava, P., & Stefanova, L. (2013). The impact of climate change on rainfall intensity-duration-frequency (IDF) curves in Alabama. *Regional Environmental Change*, 13(1), 25–33. https://doi.org/10.1007/s10113-012-0375-5

Molnar, P., & Burlando, P. (2008). Variability in the scale properties of high-resolution precipitation data in the Alpine climate of Switzerland. *Water Resources Research*, 44(10), W10404. https://doi.org/10.1029/2007WR006142

Muggeo, V. M. (2003). Estimating regression models with unknown break-points. *Statistics in Medicine*, 22(19), 3055–3071. https://doi.org/10.1002/sim.1545

Nadarajah, S., Anderson, C., & Tawn, J. (1998). Ordered multivariate extremes. *Journal of the Royal Statistical Society: Series B*, 60(2), 473–496. https://doi.org/10.1111/1467-9868.00136

Naveau, P., Huser, R., Ribereau, P., & Hannart, A. (2016). Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection. *Water Resources Research*, 52(4), 2753–2769. https://doi.org/10.1002/2015WR018552

Nhat, L. M., Tachikawa, Y., Sayama, T., & Takara, K. (2008). Estimation of sub-hourly and hourly IDF curves using scaling properties of rainfall at gauged site in Asian Pacific region. *Annual of Disaster Prevention Research Institute*. B, 51(B), 63–73. https://doi.org/10.1016/j.atmosres.2004.10.025

Ouali, D., & Cannon, A. (2018). Estimation of rainfall intensity-duration-frequency curves at ungauged locations using quantile regression methods. *Stochastic Environmental Research and Risk Assessment*, 32(10), 2821–2836. https://doi.org/10.1007/s00477-018-1564-7

Ouarda, T. B., Yousef, L. A., & Charron, C. (2019). Non-stationary intensity-duration-frequency curves integrating information concerning teleconnections and climate change. *International Journal of Climatology*, 39(4), 2306–2323. https://doi.org/10.1002/joc.5953

Over, T. M. (1995). Modeling space-time rainfall at the mesoscale using random cascades. PhD Thesis. University of Colorado.

Overeem, A., Buishand, A., & Holleman, I. (2008). Rainfall depth-duration-frequency curves and their uncertainties. *Journal of Hydrology*, 348(1–2), 124–134. https://doi.org/10.1016/j.jhydrol.2007.09.044

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

Panziera, L., Gabella, M., Germann, U., & Martius, O. (2018). A 12-year radar-based climatology of daily and sub-daily extreme precipitation over the Swiss alps. *International Journal of Climatology*, 38(10), 3749–3769. https://doi.org/10.1002/joc.5528

Paschalis, A. (2013). *Modelling the space-time structure of precipitation and its impact on basin response*. Doctoral dissertation. ETH Zurich. https://doi.org/10.3929/ETHZ-A-009917135

Ragno, E., AghaKouchak, A., Love, C. A., Cheng, L., Vahedifard, F., & Lima, C. H. (2018). Quantifying changes in future intensity-duration-frequency curves using multimodel ensemble simulations. *Water Resources Research*, 54(3), 1751–1764. https://doi.org/10.1002/2017wr021975

Renard, B., Kochanek, K., Lang, M., Garavaglia, F., Paquet, E., Neppel, L., et al. (2013). Data-based comparison of frequency analysis methods: A general framework. *Water Resources Research*, 49(2), 825–843. https://doi.org/10.1002/wcr.20087

Ritschel, C., Ulbrich, U., Névir, P., & Rust, H. W. (2017). Precipitation extremes on multiple timescales–bartlett–lewis rectangular pulse model and intensity-duration-frequency curves. *Hydrology and Earth System Sciences*, 21(12), 6501–6517. https://doi.org/10.5194/hess-21-6501-2017

Rivoire, P., Martius, O., & Naveau, P. (2021). A comparison of moderate and extreme ERA-5 daily precipitation with two observational data sets. *Earth and Space Science*, 8(4), e2020EA001633. https://doi.org/10.1029/2020EA001633

Roksvåg, T., Lutz, J., Grinde, L., Dyrrdal, A. V., & Thorarinsdottir, T. L. (2021). Consistent intensity-duration-frequency curves by post-processing of estimated Bayesian posterior quantiles. *Journal of Hydrology*, 603, 127000. https://doi.org/10.1016/j.jhydrol.2021.127000

Sane, Y., Panthou, G., Bodian, A., Vischel, T., Lebel, T., Dacosta, H., et al. (2018). Intensity-duration-frequency (IDF) rainfall curves in Senegal. *Natural Hazards and Earth System Sciences*, 18(7), 1849–1866. https://doi.org/10.5194/nhess-18-1849-2018

Savina, M., Schäppi, B., Molnar, P., Burlando, P., & Sevruk, B. (2012). Comparison of a tipping-bucket and electronic weighing precipitation gage for snowfall. *Atmospheric Research*, 103, 45–51. https://doi.org/10.1016/j.atmosres.2011.06.010

Scarrot, C., & MacDonald, A. (2012). A review of extreme value threshold estimation and uncertainty quantification. *REVSTAT–Statistical Journal*, 10(1), 33–60. https://www.ine.pt/revstat/pdf/rs120102.pdf

Schertzer, D., & Lovejoy, S. (1987). Physical modeling and analysis of rain and clouds by anisotropic scaling multiplicative processes. *Journal of Geophysical Research*, 92(D8), 22. https://doi.org/10.1029/jd092id08p09693

Sebille, Q., Fougères, A.-L., & Mercadier, C. (2017). Modeling extreme rainfall A comparative study of spatial extreme value models. *Spatial Statistics*, 21, 187–208. https://doi.org/10.1016/j.spasta.2017.06.009

Sevruk, B. (1997). Regional dependency of precipitation-altitude relationship in the Swiss Alps. In *Climatic change at high elevation sites* (pp. 123–137). Springer.

Sevruk, B., Matokova-Sadlonova, K., & Toskano, L. (1998). Topography effects on small-scale precipitation variability in the Swiss pre-Alps. *IAHS Publications-Series of Proceedings and Reports-Intern Assoc Hydrological Sciences*, 248, 8.

Sherman, C. W. (1931). Frequency and intensity of excessive rainfalls at Boston, Massachusetts. *Transactions of the American Society of Civil Engineers*, 95(1), 951–960. https://doi.org/10.1061/taceat.0004286

Sodemann, H., & Zubler, E. (2009). Seasonal and inter-annual variability of the moisture sources for Alpine precipitation during 1995-2002. *International Journal of Climatology*, 30(7), 947–961. https://doi.org/10.1002/joc.1932

Tencalicie, P., Favre, A., Naveau, P., Prieur, C., & Nicolet, G. (2020). Flexible semiparametric generalized Pareto modeling of the entire range of rainfall amount. *Environmetrics*, 31(2), e2582. https://doi.org/10.1002/env.2582

Tyralis, H., & Langousis, A. (2019). Estimation of intensity-duration-frequency curves using max-stable processes. *Stochastic Environmental Research and Risk Assessment*, 33(1), 239–252. https://doi.org/10.1007/s00477-018-1577-2

Ulrich, J., Fauer, F. S., & Rust, H. W. (2021). Modeling seasonal variations of extreme rainfall on different timescales in Germany. *Hydrology and Earth System Sciences*, 25(12), 6133–6149. https://doi.org/10.5194/hess-25-6133-2021

Ulrich, J., Jurado, O. E., Peter, M., Scheibel, M., & Rust, H. W. (2020). Estimating IDF curves consistently over durations with spatial covariates. *Water*, 12(11), 3119. https://doi.org/10.3390/w12113119
Van de Vyver, H. (2018). A multiscaling-based intensity-duration-frequency model for extreme precipitation. *Hydrological Processes*, *32*(11), 1635–1647. https://doi.org/10.1002/hyp.11516

Van de Vyver, H., & Demarée, G. R. (2010). Construction of Intensity-Duration-Frequency (IDF) curves for precipitation at Lubumbashi, Congo, under the hypothesis of inadequate data. *Hydrological Sciences Journal*, *55*(4), 555–564. https://doi.org/10.1080/02626661003747390

Veneziano, D., & Lepore, C. (2012). The scaling of temporal rainfall. *Water Resources Research*, *48*(8), W08516. https://doi.org/10.1029/2012WR012105

Veneziano, D., Lepore, C., Langousis, A., & Furcolo, P. (2007). Marginal methods of intensity-duration-frequency estimation in scaling and nonscaling rainfall. *Water Resources Research*, *43*(10), W10418. https://doi.org/10.1029/2007WR006040

Viviroli, D., Sikorska-Senoner, A. E., Evin, G., Staudinger, M., Kauzlaric, M., Chardon, J., et al. (2022). Comprehensive space-time hydrometeorological simulations for estimating very rare floods at multiple sites in a large river basin. *Natural Hazards and Earth System Sciences*, *22*(9), 2891–2920. https://doi.org/10.5194/nhess-22-2891-2022

Willems, P. (2000). Compound intensity/duration/frequency-relationships of extreme precipitation for two seasons and two storm types. *Journal of Hydrology*, *233*(1), 189–205. https://doi.org/10.1016/S0022-1694(00)00233-X

Youngman, B. D. (2019). Generalized additive models for exceedances of high thresholds with an application to return level estimation for U.S. wind gusts. *Journal of the American Statistical Association*, *114*(528), 1865–1879. https://doi.org/10.1080/01621459.2018.1529596

Youngman, B. D. (2020). evgam: An R package for generalized additive extreme value models. arXiv:2003.04067 [stat]. Retrieved from http://arxiv.org/abs/2003.04067

Yu, P.-S., Yang, T.-C., & Lin, C.-S. (2004). Regional rainfall intensity formulas based on scaling property of rainfall. *Journal of Hydrology*, *295*(1), 108–123. https://doi.org/10.1016/j.jhydrol.2004.03.003